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t — =pcoefficiente  to  the  sources  of  projectile  dispersion.  The  error 
budgets  statistically  define  projectile  dispersion  at  a level 
conducive  to  interpretation  and  comprehensive  understanding. 
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HITS  evaluates  projectile  dispersion  by  either  analytic 
Monte  Carlo  methods  employing  trajectory  equations.  The  trajectory 
equations  are  closed  form  approximations  to  the  six  dagfoes  of 
freedom  equations  of  motion.  The  closed  form  equations  presume 
the  projectile  does  not  experience  transonitT  flow  conditions 
during  any  phase  of  flight.  Thi*  limits  the  present  HITS  code 
to  this  type  of  projectile  (referred  to  here  as  "hypervelocity", 
with  no  other  connotation  intended) . The  trajectory  equations  have 
been  subjected  "to  extensive  testing.  In  all  cases  they  were  ac- 
curate t<* 'within  engineering  tolerances  . These  equations  have 
potential  application  to  operational  fire  control  computers,  the 
commutation  of  firing  tables,  and  any  other  situation  requiring 
rapid,  accurate,  low-cost  trajectory  determination. 


Since  this  report  is  a User ’3  Manual,  the  emphasis  is  placed 
on  providing  the  information  necessary  to  operate  the  code  and 
interpret  the  results.  Detailed  input/output  information  is 
summarized  in  tables  featuring  step-by-step  cook-book  instructions. 
Example  problems  are  discusser.  Appendices  present  the  theo- 
retical and  programming  fine  points,'  as  well  as  the  complete  pro- 
gram listing. 
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This  report  is  a User's  Manual  for  the  HITS  computer  code. 

The  HITS  code  statistically  evaluates  "projectile  dispersion," 
which  is  the  dispersion  associated  with  the  inflight  behavior 
of  a gun  launched  projectile.  Projectile  dispersion  is  a funda- 
mental limit  on  overall  weapon  system  effectiveness.  The  code 
is  an  automated  computer-based  analytic  procedure  for  computing 
crossrange  and  downrange  dispersion  error  budgets  and  sensitivity 
coefficients  to  the  sources  of  projectile  dispersion.  The  error 
budgets  statistically  define  projectile  dispersion  at  a level 
conducive  to  interpretation  and  comprehensive  understanding. 

HITS  evaluates  projectile  dispersion  by  either  analytic  or 
Monte  Carlo  methods  employing  trajectory  equations.  The  trajectory 
equations  are  closed  form  approximations  to  the  six  degrees  of 
freedom  equations  of  motion.  The  closed  form  equations  presume 
the  projectile  does  not  experience  transonic  flow  conditions 
during  any  phase  of  flight.  This  limits  the  present  HITS  code 
to  this  type  of  projectile  (referred  to  here  as  "hypervelocity" , 
with  no  other  connotation  intended) . The  trajectory  equations  have 
been  subjected  to  extensive  testing.  In  all  cases  they  were  ac- 
curate to  within  engineering  tolerances.  These  equations  have 


potential  application  to  operational  fire  control  computers,  the  \ 
computation  of  firing  tables,  and  any  other  situation  requiring  ] 
rapid,  accurate,  low-cost  trajectory  determination.  > 
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Since  this  report  is  a User's  Manual,  the  emphasis  is  placed  \ 
on  providing  the  information  necessary  to  operate  the  code  and  I 
interpret  the  results.  Detailed  input/output  information  is  j 
summarized  in  tables  featuring  step-by-step  cook-book  instructions. 

Example  problems  are  discussed.  Appendices  present  the  theo-  ■ 
retical  and  programming  fine  poini.s,  as  well  as  the  complete  pro-  j 
gram  listing.  i 
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1.0  ENGINEERING  PERSPECTIVE 


The  Hypervelocity  Inflight  Trajectory  Scatter  (HITS)  com- 
puter code  is  an  automated  procedure  for  evaluating  "projectile 
disperbion" , which  is  the  dispersion  associated  with  the  in- 
flight behavior  of  a gun  launched  projectile.  This  report  docu- 
ments the  code  in  the  format  of  a User's  Manual.,  Primary  emphasis 
is  placed  on  the  information  necessary  to  operate  the  code  and 
interpret  the  results.  Theoretical  aspects  and  computational  de- 
tails are  treated  in  the  appendices. 

This  chapter  addresses  the  following.  Section  1.1  identifies 
the  sources  of  projectile  dispersion  .ind,  thereby,  more  clearly 
defines  the  term.  Section  1.2  discusses  the  relationship  of  pro- 
jectile dispersion  to  overall  weapon  system  effectiveness.  ’The 
dispersion  analysis  methodology  suggested  by  the  code  is  illus- 
trated by  example  in  Section  1.3.  Section  1.4  summarizes  the 
limitations  of  the  code.  An  overview  of  the  report  is  presented 
in  Section  1.5.  The  overall  objective  of  this  chapter  is  to  en- 
able the  engineer  to  qvdckly  determine  the  utility  of  HITS  with 
respect  to  his  particular  problem, 

1. 1 Projectile  Dispersion  Sources 

The  HITS  computer  code  quantifies  projectile  dispersion. 
Projectile  dispersion  is  defined  here  to  be  the  dispersion 
associated  with  inflight  behavior  which  cannot  be  compensated 
by  fire  control.  (HITS  does  not  compute  the  dispersion  attri- 
butable to  fire  control.  However,  it  does  include  in  the  pro- 
jectile dispersion  calculation  the  effects  of  fire  control 
corrections  for  inflight  phenomena.)  Thus,  the  sources  of 
projectile  dispersion  enumerated  in  this  section  are  those 
factors  which  influence  inflight  behavior  and  cause  the  fire 
control  predicted  trajectory  to  differ  from  the  true  flight 
path. 


Table  1-1  lists  the  sources  of  projectile  dispersion.  The 
left  hand  column  contains  the  physical  mechanisms  through  which 
the  contributing  factors  of  the  right  hand  column  act.  There 
are  four  general  categories:  (1)  initial  conditions  at  the 

muzzle,  (2)  projectile  mass  properties,  (3)  projectile  aero- 
dynamic characteristics,  and  (4)  atmospheric  effects.  The 
sources  of  Table  1-1  cause  dispersion  only  insofar  as  the  fire 
control  computer  cannot  anticipate  their  effects  and  apply  cor- 
rective action.  For  instance,  winds  are  a source  of  projectile 
dispersion  if  the  weapon  system  does  not  measure  them.  If  it 
» _ , ■ 

1 "Hyper velocity"  is  used  throughout  this  report  to  indicate 
the  limitation  of  the  present  HITS  code  to  projectiles  whose 
velocity  does  not  become  transonic  at  any  point  along  the 
trajectory.  No  other  connotation  is  intended. 
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ERROR  SOURCE 


Initial  Conditions: 

Velocity  Vector 
Angle  of  Attack 
Angular  Rates 


Inertial  Characteristics: 


Weight 

Moments  of  Inertia 
Physical  Dimensions 
C.  G.  Asymmetries 


Aerodynamic  Characteristics: 


Static  Coefficients 
Dynamic  Coefficients 
Spin  Rat'? 

Trim  Angle  of  Attack 
Static  Margin 


Atmospheric  Effects: 


Winds 

Density 


TYPICAL 

CONTRIBUTING  FACTORS 


Sabot  Separation 
In-Barrel  Dynamics 
Blast  Effects 
Barrel  /ibrations 
Tipoff  Rates 
Gun  Slewing  Rates 


Manufacturing  Tolerances 


Manufacturing  Tolerances 
Ablations  Effect  on  the 
Ballistic  Coefficient 
Coefficient  Measurement 
Errors 


Temporal  and  Spatial 
Atmospheric 
Variations 


does,  only  the  wind  measurement  error  causes  projectile  dispersion. 
Although  not  mentioned  in  Table  1-1,  simplifying  assumptions  in- 
corporated in  the  fire  control  trajectory  to  minimise  computational 
requirements  are  also  considered  sources  of  projectile  dispersion. 
The  HITS  code  accounts  for  each  of  these  factors. 

1.2  Weapon  System  Effectiveness  Implications 

This  section  discusses  the  role  projectile  dispersion  plays 
in  determining  weapon  system  effectiveness  with  implications  to 
projectile  design  procedures.  The  discussion  opens  with  a brief 
overview  of  the  factors  which  affect  weapon  system  effectiveness. 

Figure  1-1  is  a block  diagram  of  a modern  gun  weapon  system. 

It  indicates  the  various  factors  that  influence  weapon  system 
effectiveness,  as  defined  by  the  probability  of  kill.  Referring 
to  Figure  1-1,  the  engagement  scenario  shown  on  the  left  hand  side 
establishes  the  geometry  of  the  encounter.  Target  motion  dynamics 
limit  target  motion  to  physically  possible  rates  and  accelerations. 
The  acquisition  sensor  continually  measures  the  target  coordinates 
and  passes  the  information  to  the  fire  control  computer.  Fire 
control's  basic  jobs  are  to  solve  the  intercept  geometry,  align 
the  gun,  and  signal  the  gunner  at  the  appropriate  time.  Its  cal- 
culations take  into  account  net  only  the  target  coordinates  but 
also  information  concerning  the  inflight  behavior  of  the  projec- 
tile, winds,  atmospheric  density,  gun  orientation,  and  various 
system  models  describing  target  dynamic  constraints,  gun  mount 
rate  and  orientation  limits,  etc.  Fire  control  issues  commands 
to  the  gun  mount  servo  system  to  point  the  gun.  The  servo  holds 
the  gun  on  the  target  while  awaiting  the  command  to  fire.  After 
the  round  has  been  fired,  the  projectile  trajectory  relative  to 
the  true  target  position  determines  the  point  of  closest  approach, 
that  is  the  miss  distance.  If  a hit  is  scored,  the  location  of  the 
hit  and  the  projectile  terminal  ballistics  determine  the  probability 
of  kill. 


Figure  1-1  Weapon  System  Effectiveness 


Projectile  dispersion  affects  overall  weapon  system  ef- 
fectiveness. The  accuracy  with  which  the  fire  control  computer 
can  predict  the  true  flight  path  of  the  projectile  poses  a 
fundamental  limitation  on  weapon  system  effectiveness.  Even 
if  the  exact  position  of  the  target  relative  to  the  gun  were 
known  and  the  target  could  be  tracked  with  infinite  precision, 
the  inability  of  fire  control  to  predict  the  flight  path  exactly 
would  cause  the  projectile  to  miss  the  target.  Thus,  projectile 
dispersion  is  a fundamental  limit  on  weapon  system  effectiveness 
because  it  defines  a level  of  accuracy  which  cannot  be  improved 
upon  by  refinements  of  other  elements  of  the  weapon  system. 

The  HITS  code  quantitatively  defines  this  limit. 

HITS  may  also  be  used  as  a projectile  design  tool.  Since  all 
projectiles  exhibit  dispersion  to  a greater  or  lesser  extent  and 
this  can  limit  overall  weapon  system  effectiveness,  it  would 
be  prudent  to  include  the  evaluation  of  projectile  dispersion 
in  all  projectile  desijn  studies.  The  HITS  code  is  an  analytic 
tool  for  performing  this  evaluation  inexpensively.  Furthermore, 
HITS  can  compute  a single  trajectory  under  a given  set  of  con- 
ditions or  under  a cyclic  permutation  of  parameters.  The  cal- 
culation is  quick,  accurate,  and  inexpensive.  Thus,  HITS  brings 
to  design  activities,  additional  capabilities  that  go  beyond  the 
assessment  of  projectile  dispersion. 

Although  projectile  dispersion  poses  a fundamental  limit 
on  weapon  system  effectiveness,  it  is  not  always  clear  what  im- 
pact projectile  dispersion  will  have  on  overall  weapon  system 
effectiveness.  Returning  to  Figure  1-1,  it  is  clear  that  over- 
all weapon  system  effectiveness  ultimately  depends  on  how  well 
the  gun  system  operates  as  a whole.  There  is  a synergistic 
effect  which  makes  the  system  better  than  the  sum  of  its  parts, 
since  one  subsystem  can  be  designed  to  compensate  for  the  errors 
of  another.  For  instance,  fire  control  could  be  designed  to 
monitor  and  correct  for  gun  pointing  errors,  as  indicated  by 
the  feedback  loop  shown  in  Figure  1-1.  Thus,  fire  control  can 
point  the  gun  better  than  an  analysis  of  the  gun-mount  servo 
system  would  indicate  was  possible.  For  this  reason,  it  is  po- 
tentially misleading  to  forecast  overall  weapon  system  effective- 
ness from  an  analysis  of  any  single  subsystem  (e.g.,  a projectile 
dispersion  analysis).  Thus,  judgment  must  be  used  in  interpreting 
projectile  dispersion  assessments  produced  by  the  HITS  code  in 
terms  of  overall  weapon  system  effectiveness. 
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1.3  Dispersion  Analysis  Methodology 


The  HITS  computer  code  suggests  a methodology  for  deter- 
mining the  susceptibility  of  projectile  designs  to  projectile 
dispersion.  This  section  demonstrates  the  two  step  method  with 
an  illustrative  example.  The  following  paragraphs  illustrate 
the  information  required,  the  role  of  the  HITS  code,  and  the 
format  of  the  results. 

Tabla  1-1  lists  the  sources  of  projectile  dispersion.  By 
analyzing  the  contributing  factors,  statistics  may  be  assigned 
to  each  source  which  quantify  the  uncertainty  as  illustrated  in 
Table  1-2:  the  error  source  model.  Displayed  in  the  left  hand 

column  of  Table  1-2  are  the  sources  of  projectile  dispersion. 

The  third  column  contains  the  nominal  values  for  the  projectile 
properties.  The  second  column  defines  the  distribution  of 
variations  about  the  nominal.  The  uncertainties  are  given  ig 
the  fourth  column:  the  standard  deviations  or  one-sigma  ( l-«r). 

values. ^ The  standard  deviations  quantify  the  level  of  uncertainty 
in  the  information  supplied  to  fire  control.  The  error  source  model 
may  be  even  more  detailed.  The  basic  error  source  model  of 
Table  1-2  could  be  expanded  to  include  correlations  between  the 
elements  of  the  error  source  model,  such  as  projectile  weight 
and  muzzle  velocity.  Whatever  the  level  of  detail,  construction 
of  the  error  source  model  concludes  the  first  step  in  evaluating 
projectile  dispersion  via  the  HITS  code. 

The  second  and  last  step  is  to  use  the  HITS  code  to  evaluate 
the  projectile  dispersion  corresponding  to  the  error  source 
model  as  suggested  by  the  schematic  block  diagram  of  Figure  1-2. 

The  code  reads  the  error  source  model  as  data  and  proceeds  to 
quantify  projectile  dispersion  in  the  form  of  an  error  budget. 

The  code  contains  two  sets  of  closed  form  trajectory  equations: 
one  represents  the  true  trajectory  and  the  other  represents  fire 
control.  A Statistical  Processor  manipulates  the  trajectory  . 
equations  in  either  an  analytical  or  Monte  Carlo  fashion  to  de- 
termine the  dispersion  statistics.  The  results  of  the  calcu- 
lations are  most  conveniently  summarized  in  the  format  of  a 
projectile  dispersion  error  budget  as  indicated  in  Figure  1-2. 

Any  simplifying  assumptions  incorporated  in  an  actual  fire  con- 
trol trajectory  are  input  to  the  code  and  taken  into  account. 

The  HITS  code  can  construct  error  budgets  in  both  cross- 
range and  downrange  directions  as  well  as  evaluate  three  di- 
mensional dispersion  indices  such  as  the  Spherical  Error  Probable 
(SEP).  Error  budgets  may  be  constructed  at  either  nominal  time  or 
nominal  range.  In  this  respect  HITS  is  very  flexible  and  can 
compute  most  every  statistic  customarily  used  to  describe  dispersion. 


1 Rayleigh/Uniform  uncertainties  are  specified  by  the  mean 
magnitude  given  in  the  third  column. 
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Table  1-2  Error  Source  Model 


ERROR  SOURCE 

STATISTICAL 

DISTRIBUTION 

MEAN 

VALUE 

STANDARD 

DEVIATION 

• INITIAL  CONDITIONS 

• VELOCITY  VECTOR 
- MAGNITUDE 

Gaussian 

11, GUO  ft/sec 

- - ....  — ■ 
1/3% 

- ORIENTATION 

Ray leigh/Un i form* ** 

3,0005  deg 

— 

• INERTIAL  ORIENTATION 
- ATTITUDE 

Rayleigh/Uniform* 

C.l  deg 

- ATTITUDE  RATE 

Rayleigh/Uniform* 

55  rad/sec 

— 

• PHYSICAL  CHARACTERISTICS 
• WEIGHT 

Gaus.ian 

0.11  lbs 

1.0% 

• MOMENTS  OF  INERTIA 
- AXIAL 

Gaussian 

1.234  x 10-6  slug-ft2 

1 2/3% 

- PITCH 

Gaussian 

2.110  x 10-5  slug-ft2 

1 2/3% 

• PHYSICAL  DIMENSIONS 
- REFERENCE  AREA 

Uniform 

3.068  x 10~3  ft2 

2/3% 

- LENGTH 

Uniform 

0.31  ft 

1/3% 

- RASE  DIAMETER 

Uniform 

0.0625  ft 

1/3% 

• AERODYNAMIC  CHARACTERISTICS 

• STATIC  COEFFICIENTS 
- DRAG  VARIATION 
EFFECTS 
VELOCITY 

Cau.'aian 

+ * 

1.0% 

ANGLE  OF  ATTACX 

Gaussian 

x.O  1/deg* 

1.0% 

- NORMAL  FOROI 

Gaussian 

1.9767  1/tad 

2.0% 

• DYNAMIC  COEFFICIENTS 
- PI-CH  DAMPING 

Gaussian 

-7.5  (none) 

20.0% 

- MAGNUS  MOMENT 

Gaussian 

0.0  (none) 

0.0 

• SPIN  RATE 

Gau..i,i<?n 

400  rad/^ec 

1.0% 

• TRIM  ANGLE  OF  ATTACK 

Rayl tigo/Uniform* 

0.1  deg 

— 

• STATIC  MARGIN 

Jdussian 

6.2a  length 

. 

1.0%  of  length 

• ATMOSPHERIC  EFFECT.' 
• CONSTANT  WINDS 

Rayleigh/Uni form* 

11.0  ft/rec 

• DENSITY  VARIATIONS 

Gaussian 

2 378  x 10“3  slug-ft3 

1.0% 

*Denotos  a Rayleigh  distribution  of  magnitude  with  a Uniform  360°  distribution 
in  orientation. 


**Drag  variation  with  velocity  closely  approximated  by  Cj^  - 
is  fit  to  two  data  points  (cx^  ' vi^  “ (0.03585,  16,740)  and 
V2)  * (0,11951,  3906).  Uncertainty  in  and 


which 


OUTPUT 


PROJF.CTI  LE 
DISPERSION 
ERROR  BUDGET 


Table  1-3  presents  the  error  budget  generated  by  HITS  for 
crossrange  dispersion  at  nominal  time  corresponding  to  the  error 
source  model  of  Table  1-2.  The  error  budget  presents  both  the 
dispersion  and  the  sensitivity  coefficient  associated  with  each 
error  source.  This  identifies  the  larger  contributors  and  those 
with  significant  potential.  The  total  dispersion  is  given  in 
the  lower  right  hand  corner.  Table  1-3  illustrates  the  level 
of  detail  and  types  of  information  available  from  HITS.  This 
is  precisely  the  level  of  detail  the  engineer  needs  to  assess 
the  situation  and  make  decisions.  Supplemental  information 
concerning  the  dispersion  probability  distribution  can  be  de- 
veloped if  desired. 

In  summary,  the  methodology  suggested  by  the  HITS  code 
consists  of  (1)  constructing  a comprehensive  error  source  model 
and  (2)  processing  it  to  obtain  a detailed  error  budget.  The 
error  budget  contains  the  information  necessary  to  evaluate  a 
projectile  with  respect  to  projectile  dispersion  program  ob- 
jectives and/or  trade-off  rival  candidate  designs.  The  method 
could  also  be  employed  in  conjunction  with  ballistic  range  tests 
devoted  to  dispersion  assessment  with  the  objective  of  deriving 
additional  insight.  Whatever  the  motive,  HITS  provides  a scien- 
tific, systematic  method  of  evaluating  projectile  dispersion. 


Presentation  of  the  error  source  model  and  error  budget, 

Tables  1-2  and  1-3,  requires  comment  concerning  their  generality. 

The  error  source  model  contains  implicit  assumptions  concerning 

(1)  manufacturing  tolerances  on  the  projectile  and  barrel, 

(2)  the  amount  of  testing  employed  to  determine  the  aerodynamic 
characteristics,  and  (3)  the  presumed  accuracy  of  wind  and  den- 
sity measurements.  The  error  budget  reflects  these  assumptions 

and  further  assumes  the  fire  control  trajectory  algorithm  contains  no 
simplifying  approximations.  Thus,  the  error  source  model  and 
error  budget  presented  here  are  only  for  the  purposes  of  illus- 
tration. They  are  not  necessarily  representative  of  any  weapon 
system  in  the  inventory, under  procurement,  or  in  development. 


1.4  Code  Limitations 


This  section  states  the  limitations  of  the  HITS  code. 

Most  of  the  items  listed  below  are  not  theoretical  limitations. 
Some  are  assumptions  made  for  modeling  convenience,  while  others 
are  merely  computer  requirements.  In  any  case,  they  are  listed 


Table  1-3  Croesrange  Dispersion  at  Nominal  Time 
(Nominal  Range  ■ 10  Kft) 


r 

ERROR  SOURCE 

UNCERTAINTY 

SENSITIVITY 

(MRAD/O 

1-  9 ERROR 
(MRAD) 

• INITIAL  CONDITIONS 
• VELOCITY  VECTOR 

- MAGNITUDE 

1/3% 

0.001913 

0.0006377 

- ORIENTATION 

0.0005  deg 

13.91 

0.006957 

• INITIAL  ORIENTATION 

- ATTITUDE 

0.1  deg 

0.01551 

0.001551 

- ATTITUDE  RATE 

55.0  rad/sec 

0.02346 

1.290 

• PHYSICAL  CHARACTERISTICS 

• WEIGHT 

1.0% 

0.009808 

0.009808 

• MOMENTS  OF  INERTIA 

- AXIAL 

1 2/3% 

* 0 

-0 

- PITCH 

1 2/3% 

-0 

- 0 

• PHYSICAL  DI’IENSIONS 

- RE  FE  REN 'IE  AREA 

2/3% 

0.003439 

0.002292 

- LENGTH 

1/3% 

0.01274 

0.004248 

- BASE  DIAMETER 

1/3% 

0.00002225 

0.000C07418 

AERODYNAMIC 

CHARACTERISTICS 

• STATIC  COEFFICIENTS 

- DRAG  VARIATION 
EFFECTS 

VELOCITY 
ANGLE  OF  ATTACK 

- NORMAL  FORCE 

• DYNAMIC  COEFFICIENTS 

- PITCH  DAMPING 

- MAGNUS  MOMENT 

• SPIN  RATE 

• TRIM  ANGLE  OF  ATTACK 

• STATIC  MARGIN 


ATMOSPHERIC  EFFECTS 

• CONSTANT  WINDS 

• DENSITY  VARIATIONS 


1.0% 

1.0% 

2.0% 

20.0% 

0.0% 

1.0% 

0.1  deg 
1%  L 


11  ft/eec 

1.0% 


-0  - 0 
0.000007418  0.000007418 

0.0008093  0.001619 

0.00001131  0.0002261 


0.0007241 

0.8103 

0.2183 


0.01538 

0.003429 


RSS  TOTAL  DISPERSION 


0.0007241 

0.08103 

0.2183 


0.1692 

0.003429 


1.322 


-9- 


below  to  aid  the  engineer  in  determining  whether  HITS  is  appli- 
cable to  a particular  problem  and,  if  so,  facilitate  installa- 
tion. The  limitations  are: 

• The  projectile  drag  model  assumes  the  projectile  does 
not  experience  transonic  flow  conditions  at  any  time 
from  launch  to  impact.  A transonic  condition  would 
cause  appreciable  error  only  if  it  persisted  over  a 
significant  portion  of  the  flight  time.  This  assump- 
tion can  be  overcome  by  additional  code  development. 
However,  the  code  as  presently  configured  is  limited 
to  "hypervelocity"  projectiles. 

§ Muzzle  blast,  sabot  separation  and  in-barrel  dynamic 
effects  are  combined  and  represented  by  equivalent 
velocity  and  attitude  perturbations  at  the  muzzle. 
Ballistic  range  data  reduction  methods  typically 
combine  these  effects  in  this  manner. 

§ Wind  gusts  and  density  variations,  that  is  fluctuations 
along  the  trajectory,  are  assumed  to  be  zero.  Co- 
variance  analysis  using  available  wind  gust  spectra 
suggest  wind  gusts  are  a negligible  source  of  dis- 
persion for  typical  projectiles.^-  On  the  other  hand, 
steady  winds  are  known  to  be  a significant  source 
and  are  accounted  for  by  the  HITS  code. 

• When  operated  in  the  Monte  Carlo  mode,  HITS  cannot 
assess  the  effects  of  error  source  model  correlations. 

A program  modification  would  be  required.  Correlation 
effects  may,  however,  be  evaluated  via  the  Analytical 
Statistical  mode. 

§ The  HITS  code  is  written  in  FORTRAN  IV  and  requires 
approximately  17 5K  bytes  of  computer  memory.  The 
card  deck  is  punched  in  the  EBCDIC  format. 


1 Gustafson,  T.  G. , Crimi,  P.,  Bellaire,  R.  G.,  "Dispersion  of 
Increased  Velocity  Pro^sctiles  - Feasibility  Phase,"  Avco 
Corporation,  AVSD-0200- 75-CR,  July  1975. 
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c The  random  number  generator  (Subroutine  RANDU)  will 
work  properly  only  on  a 32  bit  word  computer,  such 
as  an  IBM-360.  The  subroutine  can  be  easily  modi- 
fied. However,  some  care  must  be  exercised  in 
performing  the  modification  to  insure  the  generated 
random  sequences  will  reproduce  those  of  the  original 
code  on  a 32  bit  word  computer.  Otherwise,  the  full 
instructional  benefits  of  the  example  problems  pre- 
sented here  will  not  be  achieved,  since  it  will  not 
be  possible  to  reproduce  the  Monte  Carlo  simulations. 

• The  code  contains  an  interface  for  automatic  computer 
plots.  Since  plotting  software  is  hardware  specific, 
this  capability  is  not  fully  developed.  The  interface 
is  thoroughly  documented. 

• The  code  was  designed  to  facilitate  real-time  inter- 
active mode  operation  from  remote  keyboard  terminals. 
This  capability  is  latent  in  the  code,  but  not  an 
operational  reality. 

Report  Overview 


Since  this  report  is  a User's  Manual  for  the  HITS  code, 
the  emphasis  is  on  providing  the  information  necessary  to  op- 
erate the  code  and  interpret  the  results.  The  text  presents 
the  day-to-day  necessities  required  to  perform  dispersion 
analyses.  The  appendices  discuss  the  theoretical  and  programming 
aspects . 

Chapter  2 presents  an  overview  of  the  HITS  code  directed  to 
the  user.  The  objective  is  to  aquaint  the  analyst  quickly  with 
the  code.  The  basic  code  structure  and  options  are  described. 

The  trajectory  equations,  which  are  the  heart  of  the  code,  are 
discussed  from  a functional  point  of  view.  Chapter  2 concludes 
with  a collection  of  topics  peripheral  to  HITS,  but  pertinent 
to  dispersion  analysis. 

Detailed  input  information  is  supplied  in  Chapter  3.  A 
quick-reference  step-by-step  format  is  employed.  Tables  are 
used  to  define  aix  variables  and  present  all  relevant  information  ; 
comments  direct  the  user  to  sections  of  this  report  containing 
related  indepth  discussions.  Chapter  3 closes  with  an  illus- 
trative example  demonstrating  the  encoding  of  the  error  source 
model  of  Table  1-2 . 


Chapter  4 presents  four  input-output  example  problems.  Since 
the  computer  printed  output  is  self-explanatory,  the  initial  im- 
pression might  be  that  these  examples  are  superfluous.  They  are 
not.  They  illustrate  the  output  format  and  provide  numerical 
check  problems.  The  accompanying  text  stresses  interpretation. 
Thus,  users  are  encouraged  to  reproduce  these  results  to  famil- 
iarize themselves  with  HITS  and  to  verify  the  code  at  their 
facility. 

A brief  summary  is  presented  in  Chapter  5.  The  ante- 
cedants  of  the  HITS  code  and  related  analyses  are  discussed. 

The  appendices  relate  the  theoretical  and  programming 
aspects  of  the  HITS  code.  Appendix  A is  a suimnary  discussion 
of  the  relevant  aspects  of  probability  theory  and  statistics. 
Appendix  B discusses  programming  fine  points.  Appendix  C mathe- 
matically develops  the  trajectory  equations.  Appendix  D con- 
tains the  complete  program  listing. 
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2.0  OVERVIEW 


The  purpose  of  this  chapter  is  to  discuss  the  salient  aspects 
of  the  HITS  code.  The  objective  is  to  present  the  user  with  a 
clear  picture  of  the  essence  and  function  of  the  code  as  a whole. 

To  this  end,  analytical  and  programming  fine  points  are  glossed 
over  in  an  effort  to  place  maximum  emphasis  on  the  big  picture. 
Detailed  "how  to  do  it"  information  is  presented  in  later  chapters. 
Theoretical  analyses  and  programming  information  are  presented 
in  the  appendices . 

Since  projectile  dispersion  is  most  meaningfully  stated  in 
a statistical  sense,  the  HITS  code  quantifies  it  in  statistical 
terms.  Thus,  familiarity  with  the  basic  concepts  of  probability 
and  statistics  is  essential  to  (1)  the  preparation  of  the  en- 
gineering formulation  of  the  HITS  input,  (2)  a thorough  under- 
standing of  the  algorithms  employed  by  HITS,  and  (3)  the  inter- 
pretation of  the  HITS  output.  Appendix  A reviews  the  most  per- 
tinent elements  of  probability  and  statistics  with  the  objective 
of  supplying  this  understanding.  HITS  can  be  operated  by  per- 
sonnel unfamiliar  with  probability  and  statistics.  Given  the 
information  presented  in  Chapter  3,  they  can  encode  the  formulated 
problem  and  obtain  the  desired  numerical  results. 

Section  2.1  discusses  the  structure  of  the  code.  Analysis 
cptions  are  enumerated  in  Section  2.2.  The  close-form  trajectory 
equations,  which  are  the  heart  of  the  HITS  code,  are  described  in 
Section  2.3.  Section  2.4  discusses  three f topics  which  are  periph- 
eral to  HITS,  but  pertinent  to  dispersion  analysis. 

2 • 1 Code  Structure 

A flow  diagram  illustrating  the  basic  structure  of  the  HITS 
code  is  shown  in  Figure  2-1.  As  suggested  by  the  figure,  HITS 
basic  function  is  to  process  an  error  source  model  to  obtain  a 
projectile  dispersion  error  budget.  This  section  describes  the 
three  computational  modules  that  perform  this  task:  (1)  the 

Input  Processor,  (2)  the  Statistical  Processor,  and  (3)  the 
Projectile  Trajectory  Module.  The  code  contains  f facilities  for 
accessing  a Time  Phased  Data  Base  for  updating  projectile  charac- 
teristics under  direction  from  a remote  interactive  computer 
terminal;  however,,  this  capability  is  not  an  operational  reality 
at  the  present  time.  Sections  2.1.1  and  2.1.2  discuss  the  function 
of  the  Input  and  Statistical  Processors  which  were  taken  from 
previously  developed  software.  Section  2.1.3  functionally  des- 
cribes the  Projectile  Trajectory  Module,  which  was  developed 
specifically  for  dispersion  assessment. 
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HITS  Flow  Diagram 


2.1*1  Input  Processor 

The  Input  Processor  reads  error  source  model  information 
from  the  input  data  stream.  Two  operations  are  performed  once 
all  the  data  has  been  read:  (1)  the  input  data  is  checked  against 

a master  list  to  determine  whether  crucial  information  is  missing, 
and  (2)  the  experimentor  is  informed  as  to  tv*>  completeness  of  the 
stated  error  model.  Missing  data  is  automatically  filled  in  with 
preset  values. 

The  input  data  must  define  each  element  of  the  error  source 
model  (be  it  aerodynamic  coefficient,  physical  dimension,  or 
whatever)  either  as  a constant  or  as  a random  variable.  Random 
variables  must  be  defined  as  to  distribution  type,  mean  or  nominal 
value,  and  standard  deviation  (i.e.,  the  1- a uncertainty). 
Gaussian,  Rayleigh  and  uniform  random  variables  are  particularly 
easy  to  input,  but  any  distribution  is  accepted.  The  input  in- 
formation is  passed  to  the  Statistical  Processor. 

2.1.2  Statistical  Processor 

The  Statistical  Processor  derives  dispersion  statistics  by 
referring  to  the  Projectile  Trajectory  Module,  as  suggested  by 
Figure  2-1.  Each  reference  is  the  analytic  equivalent  of  a bal- 
listic range  shot  fired  under  the  conditions  specified  by  the 
Trajectory  Module  inputs.  The  Statistical  Processor  manipulates 
the  Trajectory  Module  inputs  in  accord  with  the  error  source  model 
uncertainties.  The  shot-by-shot  dispersion  returned  by  the 
Trajectory  Module  is  compiled  by  the  Statistical  Processor  to 
determine  the  dispersion  error  budget  statistics. 

2.1.3  Projectile  Trajectory  Module 

The  Projectile  Trajectory  Module  consists  of  two  sets  of 
closed  form  parametric  trajectory  equations.  The  first  computes 
the  Fire  Control  predicted  trajectory  and  the  second  simulates 
the  true  flight  path  or  "Real  World"  Trajectory.  The  Trajectory 
Module  returns  the  difference  to  the  Statistical  Processor. 

The  Fire  Control  equations  are  solved  for  the  nominal  time 
required  to  reach  nominal  range  and  the  coordinates  and  velocity 
of  the  projectile  at  that  time.  This  fixes  the  aim  point  in 
time  as  well  as  space.  Typically,  the  parameters  of  the  Fire 
Control  equations  are  set  equal  to  the  error  source  model  nominal 
values,  so  that  the  aim  point  is  based  on  nominal  muzzle  velocity, 
weight,  drag,  etc.  Simplifying  assumptions  found  in  operational 
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Fire  Control  computers  can  be  simulated  by  judiciously  selecting 
these  parameters.  Dispersion  is  computed  relative  to  the  aim 
point. 

The  Real  World  trajectory  equations  are  solved  for  the 
actual  projectile  coordinates  and  velocity  at  nominal  time  and 
at  nominal  range.  The  distance  between  the  projectile  position 
at  nominal  time  and  the  aim  point  defines  dispersion  for  engage- 
ments with  rapidly  moving  targets,  where  time  of  arrival  is  im- 
portant. The  displacement  of  the  projectile  position  at  nominal 
range  from  the  aim  point  defines  dispersion  for  engagements  with 
slowly  moving  targets.  The  parameters  of  the  Real  World  equations 
change  from  shot  to  shot  and  contain  perturbations  in  accord  with 
the  error  source  model  uncertainties.  Thus,  the  dispersion  re- 
turned to  the  Statistical  Processor  is  consistent  with  the  error 
source  model  and  incorporates  the  effects  of  Fire  Control  simpli- 
fying assumptions . 

2 . 2 Analysis  Options 

The  Statistical  Processor  of  Figure  2-1  was  extracted  from 
previously  developed  software.  As  a result,  it  has  several  de- 
sirable features  unrelated  to  statistical  dispersion  assessment. 
The  Statistical  Processor  has  four  modes  of  operation:  (1)  Single 

Trajectory,  (2)  Range  Check,  (3)  Analytical  Statistical,  and 
(4)  Monte  Carlo.  The  first  two  modes  are  of  primary  interest 
when  designing  projectiles  or  analyzing  various  launch  phenomena. 
The  latter  two  modes  are  statistical  methods  for  quantifying 
dispersion.  The  following  subsections  define  the  modes  and 
suggest  usages. 

2.2.1  Single  Trajectory 

This  mode  computes  the  dispersion  associated  with  a single 
shot  at  a given  range.  The  Single  Trajectory  mode  would  be  par- 
ticularly useful  in  selecting  nominal  design  parameters  such  as 
spin  rate.  A nominal  trim  lift  would  be  input  and  the  spin  rate 
varied  to  minimize  dispersion  at  maximum  range. 

2.2.2  Range  Check 

This  mode  evaluates  the  dispersion  associated  with  a syste- 
matic variation  of  input  parameters,  thereby  providing  valuable 
design  data.  For  instance,  a Range  Check  on  nominal  range 

(5000,  6000,  , 10,000  ft)  would  generate  a complete  trajectory. 

Alternatively,  a Range  Check  on  static  margin  (4,  o,  and  6%)  and 
spin  rate  (300,  400,  and  500  rad/sec)  would  produce  simulated 
shots  with  all  possible  combinations  of  static  margin  and  spin 
rate.  Comparison  of  the  dispersions  might  suggest  the  optimal 
combination. 


2.2.3  Analytical  Statistical 


The  Analytical  Statistical  mode  evaluates  first  and  second 
order  dispersion  sensitivity  coefficients  by  varying  the  tra- 
jectory inputs  about  the  error  source  nominal.  This  mode  com- 
putes , for  example,  the  jump  angle  sensitivity  to  cross  wind 
uncertainty;  that  is,  the  milliradians  per  ft/aoc.  The  Ana- 
lytical Statistical  mode  then  uses  the  sensitivity  coefficients 
to  evaluate  the  average  dispersion  as  well  as  the  1-v  uncertainty. 
The  dispersion  probability  distribution  is  not  determined.  The 
analytic  techniques  employed  permit  rapid,  low-cost  dispersion 
assessment. 

2.2.4  Monte  Carlo 

The  Monte  Carlo  mode  simulates  ballistic  range  tests.  It 
evaluates  not  only  dispersion  statistics  (i.e.,  mean  values  and 
a'  s)  but  also  the  dispersion  distribution.  Numerous  trajec- 
tories are  computed  using  computer  generated  random  numbers  to 
simulate  error  source  modal  uncertainties.  Each  simulated  tra- 
jectory is  comparable  to  a ballistic  range  shot,  at  a fraction 
of  the  cost.  A test  sequence  is  summarized  by  statistical  indices 
and  bar  charts  showing  the  number  of  times  the  projectiles  had 
dispersions  of  say  0 to  0.25  milliradians,  0.25  to  0.50  milli- 
radians, etc.  The  Monte  Carlo  mode  gives  the  most  complete  dis- 
persion picture. 

2 . 3 Trajectory  Equations 

The  Projectile  Trajectory  Module  contains  both  Fire  Control 
and  Real  World  closed  form  trajectory  equations.  Both  mechanize 
the  theory  presented  in  Appendix  C.  Section  2.3.1  discusses  the 
rationale  for  the  selection  of  closed  form  equations.  An  over- 
view of  the  theory  is  presented  in  Section  2.3.2.  Although  the 
trajectory  equations  are  used  here  to  define  dispersion,  there 
are  other  potential  applications,  which  include  incorporation  in 
operational  Fire  Control  computers,  the  computation  of  firing 
tables,  as  well  as  any  other  situation  requiring  rapid,  accurate, 
low-cost  trajectory  determination. 

2.3.1  Closed  Form  Rationale 

One  approach  to  determining  the  Projectile  Trajectory  would 
be  to  perform  full  numerical  simulations  of  the  six  degrees  of 
freedom  (6  DOF)  differential  equations  of  motion.  However,  since 
a Monte  Carlo  approach  requires  at  least  hundreds  and  possibly 
thousands  of  simulations  in  order  to  determine  the  dispersion 
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accurately,  the  comparatively  long  running  times  of  6 DOF  simu- 
lations makes  this  approach  impractical.  Another  approach  is  to 
develop  approximate  trajectory  equations  which  can  be  rapidly 
evaluated,  while  at  the  same  time  include  the  effects  of  the 
major  error  sources.  HITS  uses  the  latter  approach.  The  pro- 
jectile equations  of  motion  were  simplified  and  the  trajectory 
was  determined  in  closed  form.  Perturbation  equations  ate 
solved  to  refine  the  trajectory  model.  Six  degrees  of  freedom 
simulations  were  conducted  to  authenticate  the  trajectory  model 
approximations  and  are  presented  in  Appendix  C.  In  all  cases  the 
trajectory  model  checked  out  to  within  engineering  tolerances. 

2.3.2  Summary  Description 

The  trajectory  equations  are  composed  of  three  parts,  as 
illustrated  in  Figure  2-2.  The  first  part  is  a particle  tra- 
jectory which  accounts  for  a large  number  of  drag  effects  and 
the  effect  of  constant  velocity  winds.  The  second  part  is  a 
perturbation  model  which  evaluates  the  effect  of  lift  on  cross- 
ra^e  dispersion.  The  final  part  is  a perturbation  model  to 
; unt  for  the  downrange  effect  of  lift  induced  drag.  The 
pa  icle  trajectory  is  evaluated  first.  Corrections  for  lifting 
efl  ts  are  then  computed  via  the  perturbation  equations  and 
applied  to  the  particle  trajectory.  The  following  paragraphs 
discuss  the  three  components  of  the  Trajectory  Module. 


Particle  Trajectory 


The  particle  trajectory  equations  shown  in  Figure  2-2  account 
for  the  downrange  effects  of  muzzle  velocity,  projectile  weight, 
drag,  etc.,  as  well  as  crossrange  and  downrange  winds.  The 
equations  assume  a zero  angle  of  attack,  a uniform  density 
atmosphere,  and  the  absence  of  gravity,  and  a variable  drag 
coefficient.  Classically,  particle  trajectories  have  approxi- 
mated the  drag  coefficient  as  being  independent  of  velocity. 
Although  adequate  for  conventional  muzzle  velocities  and  short 
flight  times,  the  drag  coefficient  model  for  hypervelocity  pro- 
jectiles must  include  the  variation  with  velocity  in  order  to 
avoid  substantial  errors,  since  the  drag  coefficient  of  typical 
hypervelocity  projectiles  (i.e.,  slender  bodies)  can  vary  by  a 
factor  of  3 to  4 over  the  velocity  range  of  interest.  The  HITS 
particle  trajectory  models  the  drag  coefficient  as 


,1 

I 

■i 


PROJECTILE 
PROJECTILE 
INITIAL  CO 


which  can  also  be  expressed  in  terms  of  the  ballistic  coefficient 
( d = W/CdA): 


where  Cdqq  and  Kp  or  0QQ  and  K are  empirically  determined  constants. 
This  is  an  excellent  model  in  the  subsonic  and  supersonic  ve- 
locity regimes.  Ablative  effects  (i.e.,  tip  recession  and  mass 
loss)  can  be  approximated  by  making  appropriate  changes  to  the 
drag  model  parameters.  Thus,  including  the  drag  variation  with 
velocity  opens  the  door  to  considering  the  effects  of  ablation. 

Crossrange  Perturbation  Equations 

This  component  of  the  Trajectory  Model  computes  perturbations 
to  the  particle  trajectory  due  to  spurious  lift  forces  caused  by 
angle  of  attack  oscillations  and  static  trim  angles.  The  angle 
of  attack  history  is  determined.  The  effect  of  a constant  pro- 
jectile spin  rate  is  included.  Velocity  is  assumed  to  be  constant, 
so  the  crossrange  error  arising  from  transient  roll  resonance  is 
not  determined.  Aerodynamic  crossrange  perturbation  forces  are 
significant  only  during  the  first  few  pitch  oscillations.  Since 
these  occur  near  the  muzzle  while  velocity  is  essentially  constant, 
the  constant  velocity  assumption  is  realistic.  The  crossrange 
perturbations  are  applied  to  the  particle  trajectory  as  a function 
of  range.  The  particle  trajectory  includes  the  variation  in  ve- 
locity. The  solution  of  these  equations  are  applied  as  corrections 
to  the  particle  trajectory  as  indicated  in  Figure  2-2. 

Downrange  Perturbation  Equations 

The  downrange  perturbation  equations  incorporate  the  effects 
of  angle  of  attack  variations  on  the  drag  coefficient.  They  in- 
clude the  effect  of  a constant  spin  rate.  During  the  period  when 
the  projectile  is  oscillating,  the  perturbations  to  the  velocity  and 
flight  time  are  determined  by  computing  a mean  drag  coefficient 
with  and  without  the  angle  of  attack  history  as  determined  by 
the  crossrange  perturbation  equations.  The  perturbations  are 
applied  as  corrections  to  the  velocity  and  flight  time  as  deter- 
mined by  the  (zero  angle  of  attack)  particle  trajectory.  The  tra- 
jectory is  carried  beyond  the  point  of  angle  of  attack  convergence 
by  restarting  the  particle  trajectory  solution  with  the  corrected 
velocity  and  time  at  convergence  as  initial  conditions. 
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• The  radius  of  the  80%  circle  is  clone  ed  R80. 

• The  Radial  Standard  Deviation  (RSD)  is  the  rms  value 
of  the  impact  point  from  the  aim  point,  (x,  y), 

; RSD  = \f E [x2  + y2]  (2.4-2) 

where  E denotes  the  ensemble  average. 

• The  "jump  angle",  7,  is  the  average  displacement  of 
the  impact  point  from  the  aim  point.,  (x,  y), 

T * E + y2  J (2.4—3) 

when  the  impact  coordinates  x and  y are  expressed  as 
angular  deflections. 


• There  is  a 50%  probability  the  impact  point  coordinates 
relative  to  the  aim  point,  (x,  y),  will  lie  between  two 
parallel  lines  which  are  equidistant  from  the  origin 
and  are  separated  by  twice  the  Linear  Error  Probable 
(LEP). 

The  aim  point  is  usually  defined  as  the  centroid  of  the  impact 
points . 

Since  the  components  of  crossrange  dispersion  can  usually 
be  assumed  to  have  equal  magnitude  and  be  independent  Gaussian 
random  phenomenon,  the  statistical  dispersion  indices  are  pro- 
portional. The  scale  factors  are  presented  in  Figure  2-3.  The 
HITS  code  evaluates  the  crossrange  standard  deviation,  ffCR* 

The  probability  of  occurrence  associated  with  each  index  is  given 
in  Table  2-1. 

2.4.2  Correlation  Analysis  - An  Example 

The  error  source  model  of  Table  1-2  treats  each  source  as 
an  independent  random  phenomenon.  Whereas  this  is  a perfectly 
valid  assumption  in  most  instances,  the  potentiality  exists  that 
there  are  some  subtle  correlations  (i.e.,  interrelationships) 
which  could  affect  the  dispersion  statistics.  The  HITS  code  can 
evaluate  the  effects  of  correlation  provided  the  correlation  co- 
efficient can  be  determined.  This  section  illustrates  analytic 
procedures  for  evaluating  correlation  coefficients  from  design 
information.  The  most  direct  approach  would  be  to  use  ballistic 
range  test  data  to  evaluate  the  correlation.  However,  for  the 
purposes  of  illustration,  it  is  assumed  this  avenue  is  not  open 
and  the  correlation  must  be  evaluated  from  design  information. 

The  example  is  incomplete  in  the  sense  that  no  concrete  con- 
clusions are  drawn.  The  objective  is  to  illustrate  that  cor- 
relation coefficients  are  amenable  to  analyses  based  on  physical 
considerations.  The  correlation  coefficient  is  theoretically 
discussed  in  Appendix  A. 3. 3. 

Intuitively,  muzzle  velocity  variations , 5 V and  projectile 
weight  variations,  5Wp,  would  be  expected  to  be  negatively  cor- 
related to  some  extent,  since  a heavier  projectile  should  result 
in  a lower  muzzle  velocity.  These  two  sources  are  typically 
large  contributors  to  downrange  dispersion  if  they  are  considered 
uncorrelated.  However,  because  a heavier  projectile  would  be 
expected  to  slow  down  less,  it  might  be  suspected  that  when  cor- 
relation is  taken  into  account,  weight  and  muzzle  velocity  effects 
compensate,  and  result  in  signif ice ntly  less  net  downrange  dis- 
persion. Since  they  are  likely  to  be  correlated  and  their 


correlation  potentially  could  have  a substantial  effect  on  down- 
range  dispersion,  muzzle  velocity  and  projectile  weight  variations 
are  selected  for  the  purpose  of  discussing  analyses  procedures. 


Presuming  the  gun  is  a rela' ively  uniform  accelerator  suggests 
an  analysis  of  the  kinetic  energy  imparted  to  the  launch  package 
might  be  the  most  direct  approach  to  determining  the  correlation 
between  projectile  weight  and  muzzle  velocity  variations.  (An 
analysis  of  the  momentum  imparted  by  launch  might  be  a more  fruit- 
ful approach  in  a launcher  with  an  impulsive  launch  cycle.)  The 
kinetic  energy  of  the  sabot/projectile  is 
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(2.4-4) 


where  Wa  is  the  sabot  weight  and  g is  the  acceleration  due  to 
gravity.  Viewing  weight  as  the  dependent  variable,  Eq.  (2.4-4) 
can  be  linearized  about  the  nominal  launch  conditions  and  solved 
for  the  variation  in  velocity 


fiv  i \ 


<1  - •>) 


(2.4-5) 


where  f?  is  the  design  ratio  of  projectile  weight  to  the  total  package 
weight.  Equation  (2.4-5)  determines  the  muzzle  velocity-projectile 
weight  correlation  coefficient.  It  is 
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(2.4-6) 


The  correlation  coefficient  pWgWp  is  the  correlation  coefficient 
of  the  projectile  and  sabot  weights.  Since  these  are  the  result 
c . different  manufacturing  processes,  Pwswp  is  most  probably 
zero.  and  are  the  standard  deviations  of  the  projectile 

weight  and  the  muzzle  velocity  and  are  given  by  the  basic  error 
source  model.  is  the  standard  deviation  of  the  kinetic 

energy  imparted  to  the  launch  package.  Assuming  the  launch  cycle 
is  efficient,  should  be  dominated  by  the  variations  in  the 

chemical  energy  released  by  the  propellant.  (Other  factors  in- 
fluencing <7^  would  include  heat  absorbed  by  the  gun,  the  energy 
lo8 t by  in-barrel  balloting  and  friction,  and  the  residual  energy 
in  the  gun  gases.)  The  launch  package  kinetic  energy-projectile 
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weight  correlation  coefficient,  pKeWp,  is  an  interesting  one. 

At  first  these  two  quantities  would  appear  to  be  unrelated 
( pKeWp  = 0).  However,  the  projectile  weight  has  an  influence 
on  the  diameter  of  the  projectile  which  in  turn  affects  the  tight- 
ness of  the  seal  between  the  launch  package  and  the  barrel.  By 
this  argument  an  increased  weight  should  cause  more  kinetic  en- 
ergy to  be  imparted  to  the  launch  package.  Analysis  of  the  pro- 
jectile geometry  and  consideration  of  typical  variations  in 
material  mas'  densities  would  determine  the  appropriate  value  of 
pKeWp.  All  values  would  be  substituted  into  Eq.  (2.4-5)  to  de- 
termine the  projectile  weight-muzzle  velocity  correlation 
coefficient . 

This  example  has  illustrated  the  determination  of  correlation 
coefficients.  They  may  be  evaluated  from  ballistic  range  test 
data  or  design  analyses.  The  analytical  procedures  are  straight- 
forward, systematic,  and  scientific  and  are  based  on  the  first 
principles  of  the  physical  sciences . 

2.4.3  Target-Fixed  Crossrange  Dispersion 

The  HITS  code  defines  projectile  dispersion  in  a coordinate 
frame  fixed  with  respect  to  the  gun  implacement.  The  dispersion 
must  be  transformed  into  a frame  moving  with  the  target  in  order 
to  assess  weapon  system  effectiveness.  The  conversion  to  target- 
fixed  coordinates  is  usually  performed  during  an  engagement  analysis 
and  can  be  quite  elaborate.  This  section  presents  a simplified 
procedure  which  allows  the  analyst  to  quickly  det  .rmine  the  approxi- 
mate target-fixed  crossrange  dispersion. 
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For  the  intercept  geometry  of  Figure  2-4,  the  target-fixed 
crossrange  dispersion  in  milliradians , ffTCD'  is  given  by 


7TCD  = VCTCr2  + (Sc/D  aDR)2 


(2.4-7) 


where  and  are  the  crossrange  and  dowr.range  projectile 

dispersion  standard  deviations  in  milliradians  and  percent  of 
range,  respectively,  and 


where  Vp/Vt  is  the  projectile  to  target  velocity  ratio  in  the 
vicinity  of  the  target,  and  ^ is  the  deflection  angle,  i.e.,  the 
angle  between  the  target  velocity  vector  and  the  bore  sight. 

The  coefficient  S^/d  is  the  equivalent  crossrange  dispersion 
sensitivity  to  downr-nge  dispersion.  The  sensitivity  coefficient 
is  evaluated  in  Figure  2-5  for  representative  velocity  ratios. 
Maximum  sensitivity  occurs  at  a deflection  angle  of 


^«n«*  = 90°  + 


(2.4-9) 


The  maximum  sensitivity  is 


(£C/D> 


max 


10 

y(vp/vt)2  - 1 


(2.4-10) 


and  can  be  used  to  perform  a worst-case  analysis. 

This  section  has  presented  a simple  approximate  technique 
for  converting  projectile  dispersion  into  a target-fixed  coordinate 
frame.  The  results  may  be  used  in  conjunction  with  target  vulner- 
ability, projectile  terminal  ballistics,  and  target  lethality 
data  to  estimate  the  effect  of  projectile  dispersion  on  weapon 
system  effectiveness. 
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3.0  INPUT  DESCRIPTION 


This  chapter  describes  the  order  and  content  of  the  input 
cards.  The  input  for  a single  case  consists  of  three  card  groups. 
Section  3.1  defines  the  Statistical  Processor  Control  Card  (Card 
Group  1).  Input  Processor  Control  Cards  (Card  Group  2)  are  dis- 
cussed in  Section  3.2.  The  Optional  Input  Processor  Control 
Cards  (Card  Group  3)  are  described  in  Section  3.3.  The  card 
groups  are  defined  by  cook-book  tables  that  give  the  card  format 
and  options.  Comments  in  the  tables  alert  the  user  to  related, 
indepth  discussions  presented  in  other  portions  of  the  text. 
Section  3.4  illustrates  encoding  of  an  error  source  model. 

The  card  groups  for  a single  case  must  appear  in  the  order 
suggested  by  the  card  group  numbers.  The  card  formats  utilize 
fixed  numeric  fields  to  facilitate  da^n  base  operations.  IBM- 
029  keypunch  control  cards  are  presen  ad  to  simplify  keypunching. 
All  integer  variables  must  be  right-hand  justified.  Multiple 
cases  may  be  stacked.  Each  case  in  the  stack  must  be  complete 
in  itself  and  not  simply  changes  from  the  prior  case. 

3 . 1 Statistical  Processor  Control  Card  (Card  Group  1) 

The  Statistical  Processor  Control  Card  is  the  first  card 
appearing  in  the  deck  describing  a single  case.  This  card 
directs  the  Statistical  Processor  to  enter  a specific  mode  of 
operation.  Table  3-1  defines  the  format  of  the  card  and  the 
various  options.  Figure  3-1  illustrates  the  associated  IBM-029 
keypunch  control  card. 

The  Statistical  Processor  Control  Card  (Card  Group  1) 
described  in  Table  3-1,  defines  13  variables.  Only  I0PRNT  and 
MC0PT  play  significant  roles  and  are  discussed  in  depth.  I0PRNT 
and  MC0PT  control  thei  mode  of  the  Statistical  Processor  as  de- 
scribed in  the  following  two  paragraphs. 

Referring  to  Table  3-1,  the  I0PRNT  variable  serves  a dual 
purpose.  When  I0PRNT  = 0,  the  Statistical  Processor  identifies 
the  case  as  either  a Single  Trajectory  or  Range  Check  simulation, 
and  relinquishes  control  to  the  Input  Processor.  The  Input 
Processor  determines  the  proper  mode  from  the  contents  of  Card 
Group  2.  Whenever  I0PRNT  > 0,  the  Statistical  Processor  enters 
the  Analytical  Statistical  mode  and  requests  the  Input  Processor 
to  supply  the  contents  of  Card  Group  2. 
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The  MC0PT  variable  of  Table  3-1  determines  whether  or  not 
the  Statistical  Processor  will  enter  the  Monte  Carlo  mode  after 
completing  Analytical  Statistical  mode  calculations.  If  MC0PT 
» 1,  the  Statistical  Processor  automatically  shifts  into  the 
Monte  Carlo  mode  after  the  Analytical  Statistical  mode  results 
have  been  used  to  set  up  the  histograms. 

3 . 2 Input  Processor  Control  Cards  (Card  Group  2) 

The  Input  Processor  Control  Cards  appear  second  in  the  deck 
describing  a single  case.  Each  card  directs  the  Input  Processor 
to  assign  qualities  to  a trajectory  variable.  These  qualities 
affect  the  calculations  performed  as  described  in  Section  3.2.1. 
The  trajectory  variables  are  defined  in  Section  3.2.2. 

3.2.1  Card  Group  2 Assigned  Qualities 

Each  Group  2 card  defines  two  parameters:  C0DE#  and  TYPE. 

C0DE#  (i.e.,  code  number)  identifies  the  trajectory  variable  to 
which  the  qualities  determined  by  the  value  of  TYPE  are  to  be 
assigned.  This  is  one  of  the  more  interesting  aspects  of  the 
HITS  code.  There  are  two  basically  different  TYPE  specifications 
that  can  be  assigned.  They  are  described  in  the  succeeding 
paragraph.  It  is  assumed  the  user  wishes  to  determine  the  ef- 
fect of  projectile  weight  variations  on  downrange  dispersion  at 
nominal  time,  for  the  purposes  of  discussion. 

Whenever  1 < TYPE  < 5,  the  variable  is  assigned  the  dis- 
tinction of  being  an  independent  variable,  that  is,  an  independent 
variable  whose  value  is  determined  by  input  data  rather  than 
preset  (or  default)  values.  In  the  example,  projectile  weight 
is  the  independent  variable  and  would  be  assigned  a TYPE  between 
1 and  5 depending  on  whether  it  is  to  be  treated  as  a determin- 
istic constant  (different  than  the  preset),  a range  check  var- 
iable, or  a random  variable.  Whenever  7 < TYPE  < 8,  the 
variable  is  defined  to  be  a dependent  variable,  that  is,  a de- 
pendent variable  for  which  printed  output  is  to  be  produced. 

In  the  example,  downrange  dispersion  at  nominal  time  would  be 
declared  a TYPE  = 7 or  TYPE  = 8 variable,  since  the  user  desires 
to  have  the  dispersion  displayed.  Thus,  the  C0DE#  identifies 
the  variable  and  TYPE  determines  what  is  to  be  done  with  it. 

All  independent  variables  not  mentioned  in  Card  Group  2 are 
given  their  preset  values.  All  dependent  variables  not  men- 
tioned are  computed  but  not  printed. 

A Group  2 card  may  require  more  information  than  just  the 
C0DE#  and  TYPE  to  completely  define  the  variable.  Table  3-2 
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Table  3-2  Input  Processor  Control  Cards  (Card  Group  2)  (cont‘d) 


ST  is  defined  later  in  this  ttbl 


Table  3-2  Input  Processor  Control  Cards  (Card  Group  2)  (cont'd) 


blank. 


Table  3-2  Input  Processor  Control  Cards  (Card  Group  2 (concl'd) 
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TYPE  = 4 NT  is  the  number  of  pairs  of  points  • For  TYPE  = 4,  NT  is  the  number 

used  to  represent  the  probability  ordinate-abscissa  pairs. 


explains  the  details.  In  Table  3-2  hash  marks  are  used  to  de- 
limit the  extent  of  the  comments  made  in  the  right  hand  column. 
Comments  identified  by  stars  command  action  on  the  part  of  the 
user.  Group  2 cards  may  be  arranged  in  any  order.  IBM-029  key- 
punch control  cards  for  Group  2 are  presented  in  Figure  3-2. 

3.2.2  Trajectory  Variable  Code  Numbers 

HITS  uses  "code  numbers"  to  identify  variables  in  input 
data  cards  and  the  printed  output.  C0DE#  in  Card  Group  2 is 
just  the  first  occurrence  of  code  numbers  in  the  discussion. 

This  section  discusses  Table  3-3  which  catalogs  the  code  numbers. 

The  code  number  of  a variable  is  its  address  in  the  HITS' 
active  storage  array  (i.e.,  the  0E  array).  Table  3-3  lists  the 
code  numbers  and  their  associated  FORTRAN  names,  preset  values, 
units,  and  analysis  names.  The  code  numbers  also  link  the  HITS 
variables  to  the  closed  form  trajectory  equations  of  Appendix  C 
via  Table  C-l.  An  asterisk  attached  to  a FORTRAN  name  denotes  a 
potential  independent  variable.  All  others  may  be  defined  as  de- 
pendent variables.  The  presets  are  the  nominal  values  of  the 
Table  1-2  error  source  model. 

Referring  to  Table  3-3,  the  code  numbers  from  1 to  100  de- 
fine Fire  Control  trajectory  variables.  Code  numbers  from  101 
to  200  are  the  Real  World  trajectory  variables.  Variables  with 
code  numbers  from  201  to  300  are  system  controls.  These  are 
particularly  important  since  they  can  greatly  reduce  the  amount 
of  input  data  and  the  number  of  computer  runs.  Note  that  vari- 
able 203,  (i.e.,  IFC)  distinguishes  between  Real  World  and  Fire 
Control  trajectory  solutions.  Variables  301  to  400  are  measures 
of  dispersion.  Variables  401  to  500  are  computed  trajectory 
quantities.  Code  numbers  from  501  to  600  denote  variables  de- 
fining the  terminal  conditions  of  the  Real  World  and  Fire  Control 
traj ectories. 

3.3  Optional  Input  Processor  Control  Cards  (Card  Group  3) 

Optional  Input  Processor  Control  Cards  (Card  Group  3) 
appear  third  (last)  in  the  deck  describing  a single  case.  Card 
Group  3 may  not  appear  when  the  Statistical  Processor  is  operated 
in  either  the  Single  Trajectory  or  Range  Check  modes.  It  has  no 
effect  on  the  Monte  Carlo  mode.  Only  the  Analytical  Statistical 
mode  is  affected  by  Group  3 data. 

Each  Group  3 card  directs  the  Analytical  Statistical  Pro- 
cessor to  treat  two  of  the  statistically  defined  independent 
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Figure  3-2  IBM-029  Keypunch  Control  Cards  for  Card  Group 


Table  3-3  HITS  Trajectory  Parameter  List 


Code  FORTRAN 
Number  Name 


Preset 


11,000. 


Units 


ft/sec 


Remarks 


RH0F* 


CXlF* 

VlF* 

CX2F* 

V2F* 


0. 


2.378xl0"3 


3.585xl0“2 
1.674xlC4 
1.195X10"1 
3 . 906xl03 


ft/sec 

ft/sec 


F.C.  'Nozzle  Velocity 


F.C.  Steady  Winds: 
Down  Range 
Horiz.  Cross  Wind 


O 

slugs-ftJ  F.C.  Atmospheric  Density 


ft/sec 


ft/sec 


F.C.  Projectile  Drag  Model 

'a  - Cxoo  + *^v2 

Fit  to 

(CXlFf  VlF)  & CX2F, V2F) 
CX2F  < CXlF 

VX2F  > VX1F 


THETOF* 

PSIOF* 

TD0TOF* 

PD0TOF* 


deg 

deg 

rad/sec 

rad/sec 


F.C.  Initial  Attitude 
Vertical  Plane 


Horizontal  Plane 


♦Denotes  Potential  Independent  Variable 
(1) Indicates  Fire  Control  Parameter 


Table  3-3  HITS  Trajectory 


Code 

Number 

FORTRAN 

Name 

20 

GAM OF* 

21 

AZOF* 

22 

AF* 

23 

ELLF* 

24 

DF* 

25 

WF* 

26 

AIXF* 

27 

AIYF* 

28 

PF* 

29 

ALPC0F* 

30 

CAAF* 

Preset 


Units 


3 . 068x10“ 2 ft2 
3 . 1x10“ 1 ft 

6.25xl0-2  ft 
l.lxlO'1  lbs 

1.235xl0-6  slugs- ft2 
2.110xl0“5  slugs-ft2 


400. 


rad/sec 


deg 


1/deg-1 


♦Potential  Independent  Variable 


rameter  List  (Cofct'd) 


Remarks 


F.C. 

Initial  Attitude 

(Cont'd) 

Vo 

y2 

r> 

F.C. 

A 

Projectile  Physical 
Characteristics: 

- Base  Area  ( 7r  D2/4 ) 

L 

- Length 

D 

- Base  Diameter 

W 

- Weight 

Ix-  Roll  Moment  of 

Inertia 

I 

y-  Pitch  Moment  of 

Inertia 

F.C. 

Spin  Rate 

F.C. 

Angle  of  Attack 
Oscillation  Convergence 
Criterion 

F.C.  Second  Order  Drag  Effect 
ACV  = CAAF  a2 


r 


Table  3-3  HITS  Trajectory  Parameter  List  (Cbttt'd) 


Code 

Number 

FORTRAN 

Name 

Preset 

Units 

Remarks 

101 

VO* 

11,000 

ft/aec 

R.W.^Wzzle  Velocity 

m 

WX* 

wz* 

0. 

0. 

ft/sec 

ft/sec 

R.W.  Steady  Winds 
Down  Range 
Horiz.  Cross  Wind 

104 

RH0* 

2 . 378x10"" 3 

slugs-ft3 

R.W.  Atmospheric  Density 

105 

CXI* 

3.585xl0“2 

— 

R.W.  Projectile  Drag  Model 

Cx  - ^ + K/V* 

106 

107 

VI* 

CX2* 

1.674xl04 

1.195xl0_1 

ft/sec 

Fit  to 

(CXI, VI)  & (CX2,V2) 
CX2  < CXI 

108 

V2* 

3.906X103 

ft/sec 

VX2  > VX1 

109 

XN* 

10,000 

ft 

R.W.  Nominal  Range 

R.W.  Aerodynamic  Effects: 

110 

CNA* 

1.9767 

1/rad 

• °Na 

111 

SM* 

0.062 

e Static  Margin 

112 

CKQ* 

-7.5 

— 

* 

1.13 

CMPA* 

0. 

• Magnus  Moment 

114 

ATRMS* 

0. 

deg 

e Static  Trim  ((^=0) 
-Angle  of  Attack 

115 

BTRMS* 

0. 

deg 

-Angle  of  Sideslip 

R.W.  Initial  Attitude 

116 

THETO* 

0. 

deg 

Vertical  Plane 

\A 

117 

PSIO* 

0. 

deg 

*• 

X° 

V 

118 

TD0TO* 

0. 

rad/sec 

▼ 1 

0o  Horizontal  Plane  T 

^ . 1 

A 

* 

119 

P®0Tt>* 

0. 

rad/sec 
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^Potential  Independent  Variable 
(l)-Denotes  Rteal  World  Parameter 


Table  3-3  HITS  Trajectory  Parameter  List  (Cont'd) 


Code 

FORTRAN 

1 

Number 

Name 

Preset 

Units 

Remarks 

R.W.  Initial  Attitude  (Cont'd) 

120 

GAMO* 

0. 

deg 

>0 

121 

AZO* 

0. 

deg 

yz 

0 

- 

R.W.  Projectile  Physical 

Characteristics 

122 

A* 

3 . 068xl0-3 

ft2 

A - Base  Area  ( * D~/4) 

123 

ELL* 

3 . lxlO-1 

ft 

L - Length 

124 

D* 

6.25xl0-2 

ft 

D - Base  Diameter 

125 

W* 

l.lxlO"1 

lbs 

W - Weight 

126 

AIX* 

1.235xl0“6 

slugs- ft2 

Ix-  Roll  Moment  of  Inertia 

127 

AIY* 

2.110xl0~5 

slugs-ft2 

Iy-Pitch  Moment  of  Inertia 

128 

P* 

400. 

rad/sec 

R.W.  Spin  Rate 

129 

ALPC0N* 

0.5 

deg 

R.W.  Angle  of  Attack 

Oscillation  Convergence 
Criterion 

R.W.  Second  Order  Drag  Effect 

130 

CAA* 

0. 

1/deg2 

ACX  = CAA  a2 

Externally  Supplied  SEP  Center 

131 

XBART* 

0. 

ft 

^ •• 

*T  (Displacements  from 

Fire  Control  Nominal 
Aim  Point) 

132 

YBART* 

0. 

ft 

133 

ZBART* 

0. 

ft 

5T-J-  • 

Externally  Supplied  CEP  Center 

134 

YBARX* 

0. 

ft 

sy  (Displacements  from 

Fire  Control  Nominal 

135 

ZBARX* 

0. 

ft 

Aim  Point) 

♦Potential  Independent  Variable 


Table  3-3  HITS  Trajectory  Parameter  List  (Cont'd) 


Code  FORTRAN 
Number  Name 


201  I 


Preset 


Units 


Remarks 


System  Control 


IFCRW=1  F.C.  Parameters  are 
sat  equal  to  R.W. 
nominal  values 


202  IFCRC* 


IFCRW=>0  F.C.  Parameters 

determined  solely 
by  inputs 


System  Control 


IFCRC=1  F.C.  Parameters  are 
set  equal  to  R.W. 

, values  for  each 

Range  Check 
Combination 


203  IFC 


204  IDUMP* 


IFCRC=0  F.C.  Parameters 

determined  solely 
by  inputs. 


Internal  Sequencing  Variable 

IFC  =1  F.C.  Solution 
IFC  «=  0 R.W.  Solution 


System  Print  Control 


IDUMP  = 1 Print  lots  of 

intermediate  data 


205  ICONV* 


IDUMP  = 0 Print  summary  data 
only 


Trajectory  Module  Iteration 
Control 


TOL  = 10. 


ICONV 


206  ICNCL* 


* Potential  Independent  Variable 


Trajectory  Aerodynamic 
Control  Parameter 

ICNCL  = 1 Caused  in 
croscrange  perturbation. 
Recommended  for  most 
applications . 

ICNCL  *=  0 Caused  in 
crossrange  perturbation. 


Table  3-3  HITS  Trajectory  Parameter  List  (Cont'd) 


Code 

FORTRAN 

Number 

Name 

Preset 

Units 

Remarks 

Internally  Computed  SEP  Center 

301 

XBARTI 

— 

ft 

STt  (Displacements  from  Fire 

Control  Nominal  Aim  Point) 

302 

YBARTI 

— 

ft 

SJt 

303 

ZBARTI 

ft 

8?t 

Internally  Computed  CEP  Center 

304 

YBARXI 

-mmm 

ft 

Sj  (Displacements  from  Fire 
_ Control  Nominal  Aim  Point) 

30b 

ZBARXI 

— 

ft 

8zx 

Projectile  Position  at 
Nominal  Time 

306 

DXT 

— 

ft 

8*t  (Displacement  from  Fire 

Control  Nominal  Aim  Point) 

307 

DYT 

— 

ft 

SVt 

308 

DZT 

— 

ft 

8zt 

Projectile  Position  at 
Nominal  Range 

309 

DYX 

— 

ft 

1 (Displacement  from  Fire 

x Control  Nominal  Aim  Point) 

310 

DZX 

ft 

Velocity  Difference  at 
Nominal  Time 

311 

DVXT 

— 

ft/sec 

5v*t  (Difference  from  F.C. 

Value  at  Nominal  Time) 

312 

DVYT 

ft/sec 

b Vyf 

313 

DVZT 

ft/sec 

5VZt 

Velocity  Difference  at 
Nominal  Range 

314 

DVXX 

-- 

ft/sec 

SVXX  Difference  from  F.C. 

315 

DVYX 

— 

ft/sec 

. Value  at  Nominal  Time) 

0 Vyjj 

316 

DVZX 

— 

ft/sec 

«vzx 

317 

DT 

— 

sec 

Time  to  Nominal  Range  less 

F.C.  Nominal  Time  to 
Nominal  Range 

6- 


Table  3-3  HITS  Trajectory  Parameter  Diet-  (Cowfr**}  ' 


Code 

FORTRAN 

Number 

Name 

Preset 

Units 

Remarks 

318 

RADT+ 

— 

ft 

Radial  Displacement  from 
Externally  Supplied  SEP  Center 

319 



RADX+ 

— 

Radial  Displacement  from 
Externally  Supplied  CEP  Center 

320 

RH0T(1J 

— 

ft 

Radial  Displacement  from 
Internally  Computed  SEP  Center 

321 

„ (2) 

RH0X 

ft 

Radial  Displacement  from 
Internally  Computed  CEP  Center 

Projectile  Position  at 
Nominal  Time 

32; 

323 

DXTI 

ft 

Sxti  (Displacement  from  Internally 

Computed  SEP  Center) 

®yti 

DYTI 

— 

ft 

324 

DZTI 

ft 

8zti 

325 

Projectile  Position  at 
Nominal  Time 

DXTE 

ft  ! 

s*te  (Displacement  from  Externally 
R Supplied  SEP  Center) 

326 

DYTE 

ft 

327 

DZTE 

ft 

5zte 

328 

Projectile  Position  at 
Nominal  Range 

DYXI 

— 

ft  ; 

(Displacement  from  Internally 

329 

DZXI 

ft  i 

§2^.  Computed  CEP  Center) 

1 

Projectile  Position  at 
Nominal  Range 

DYXE 

— 

ft 

<5yxe  (Displacement  from  Externally 

DZXE 

— 

ft 

8zxe  Supplied  CEP  Center) 

Correlations  of  Projectile 

Position  at  Nominal  Time 

332 

DXDYT 

— 

ft2 

i 

SxtSyt  (Displacements  from 

333 

DXDZT 

— 

rt2  1 

5*t5zt  Fire  Control  Nominal 

334 

DYDZT 

— 

ft2 

8yt  S r.(  Aim  Point) 

"*  v’YPE  = 8 only 

(1) Requires  Code  Numbers  306,  307,  and  308  be  TYPE  = 7 or  TYPE  = 8. 
^2) Requires  Code  Numbers  309  and  310  be  TYPE  * 7 or  TYPE  = 8. 


Table  3-3  HITS  Trajectory  Parameter  List  (Cont'd) 


Code 

Number 

FORTRAN 

Name 

Preset 

Units 

Remarks 

335 

DYDZX 

ft2 

Correlation  at  Nominal  Range 
8yx8zx  (Displacements  from 

F.C.  Nominal  Aim 
Point) 

336 

337 

338 

DXDYTI 

DXDZTI 

DYDZTI 

— — 

ft2 

ft2 

Correlations  about  Internally 
Computed  SEP  Center  at  Nominal 
Time 

3xtisVti  (Displacements  from 

Sx  ■ 8z  Internally  Computed 

11  M SEP  Center) 

syti  ^2ti 

339 

DYDZXI 

— 

ft2 

Correlations  about  Internally 
Computed  CEP  Center  at  Nominal 
Range 

$yxi8zxi  (Displacements  from 

Internally  Computed 
CEP  Center) 

340 

341 

342 

DXDYTE 

DXDZTE 

DYDZTE 

— 

ft2 

ft2 

ft2 

Correlations  about  Externally 
Supplied  SEP  Center  at  Nominal 
Time 

Sx^Syte  (Displacements  from 

8xte8z te  Externally  Supplied 

SEP  Center) 

5yte  8zte 

343 

DYDZXE 

ft2 

Correlations  about  Externally 
Supplied  CEP  Center  at  Nominal 
Range 

5 yxe  8 zxe  (Displacements  from 
Externally  Supplied 
CEP  center) 

343 


DYDZXE 


(Displacements  from 
Externally  Supplied 
CEP  center) 


Pable  3-3  BITS  Trajectory  Parameter  List  (Oont'dj 


Code  FORTRAN 
Number  Name 


400  ALPHA 

401  at*t 

402  ALPTRM 

403  ADPIJD0 


BETTRM 


BETA0 


BETAD0 


408  BETA 


409  ALPMAX 


410  ALPMIN 


411  CAYD 


412  CD8 


CMTHTD 


CMTHTA 


CAYl 


CAY  2 


CDA0B 


419  CDAB 


Preaet 


Units 


Remarks 


a (t^)  01  ° 
a (tN)  or  a (*n^ 

Rolling  Trim  a 


rad/sec 


Rolling  Trim  0 


rad/sec 


0 (tN)  or  0(*N) 


“UPPER 


“LOWER 


(ft/seo)  Kj) 


1/rad 


Cn  (« 


Table  3-3  HITS  Trajectory  Parameter  List  (Cont'd) 


Code 

Number 

420 

421 

422 

423 

424 

425 

426 

427 

428 

429 

430 
4 31 

432 

433 
4 34 

435 

436 

437 

438 


PHIO 


PHIl 


439 


PHI  2 


:ad 


^2 


Table  3-3  HITS  Trajectory  Parameter  List  (Cont'd) 


Code 

FORTRAN 

Number 

Name 

Preset 

Units 

Remarks 

460 

WO 

-- 

sec-1 

“’o 

461 

W1 

— 

sec-'*' 

o)l 

462 

W2 

_ _ 

sec-^ 

Cl>  2 

463 

WL2 

— ■ mm 

sec-2 

' 2(o02  + AX2) 

464 

XLAMO 

mm  mm 

sec-*- 

A 

o 

465 

XLAM1 

— 

sec"*1 

Xl 

466 

XLAM2 

— 

sec"*- 

*2 

467 

XNUl 

— 

rad 

•'l 

468 

XNU2 

— 

rad 

i'2 

469 

XG 

— 

ft 

XG 

470 

XMU 

— 

— 

V- 

471 

XJAY 

rad 

JAy 

472 

XJAZ 

— 

rad 

JAz 

473 

XJA 

— 

rad 

JA 

474 

YA 

— 

ft 

*a 

475 

ZA 

— 

ft 

Za 

476 

ALPHO 

— 

rad 

ao 

tl 


Table  3-3  HITS  Trajectory  Parameter  List  (Concl'd) 


Code  FORTRAN 
Number  Name 


Preset 


Units 


Remarks 


XTFC 

YTFC 

ZTFC 


VXTFC 

VYTFC 

VZTFC 


ft/sec 

ft/sec 

ft/sec 


ft/sec 

ft/sec 

ft/sec 


F.C.  Established  Nominal  Time 
at  Nominal  Range 


F.C.  Computed  Position  at 
Nominal  Time 

*fc  ((n)  (Displacements  from 

„<  <.N)  Mu*ile> 



F.C.  Computed  Velocity  at 
Nomine  1 Time 

V ('N> 

V*> 

' 

Projectile  Position  at 
Nominal  Time 

* (Displacements  from 

y (tN)  Muzzle) 

*<«N> 

Projectile  Velocity  at 
Nominal  Time 
v,  <tN) 

vy  (tN) 

V*  <‘N> 

Projectile  Position  at 
Nominal  Range 

y (*n^  (Displacements  from 
z(xN)  Muzzle) 

Projectile  Velocity  at 
Nominal  Range 

V*  <*N> 

Vy  (ZN) 

^ w 

Projectile  Time  to  Reach 
Nominal  Rang  ^ 


variables  as  being  correlated.  Thus  each  variable  appearing  in 
Card  Group  3 must  also  appear  in  Card  Group  2 witn  a TYPE  speci- 
fication of  2,  3,  or  4.  Table  3-4  defines  the  format  of  Group  3 
cards  and  provides  pertinent  details.  The  comments  denoted  by 
stars  command  action  on  the  part  of  the  user.  The  IBM-029  key- 
punch control  cards  are  presented  in  Figure  3-3.  Group  3 cards 
may  appear  in  any  order. 

3.4  Error  Source  Model  Encoding 

This  section  illustrates  the  encoding  of  the  error  source 
model  of  Table  1-2 . Primary  attention  is  focused  on  Card  Group 
2.  Maximum  instructional  benefit  will  be  achieved  by  the  user 
who  works  the  problem  independently  and  refers  to  the  text  simply 
for  verification. 

The  first  task  is  to  identify  the  code  numbers  (i.e., 
trajectory  variable  names),  corresponding  to  line  items  in  the 
error  source  model.  Table  3-5  is  a complete  listing  arrived  at 
by  cross  referencing  the  error  source  model.  Table  1-2,  and  the 
Parameter  list.  Table  3-3.  These  are  Real  World  trajectory 
variables.  Each  of  the  Rayleigh/Uniform  distribut  d sources 
are  two-dimensional  and,  as  a result,  are  describee  by  two  tra- 
jectory variables. 

With  the  code  numbers  in  hand.  Card  Group  2 may  be  con- 
structed by  reference  to  Table  3-2.  For  instance. 

• Muzzle  velocity  ( C0DE#  = 101)  is  a Gaussian 
distributed  uncertainty  (TYPE  = 2).  The 
nominal  or  mean  value  is 

4 

VALUE  = 1.1  x 10  ft/sec 
The  standard  deviation  is 

STDEV  = (1/3%) (1.1  x 104)  = 3.666  x 101  ft/sec 
and 

T0L  = 3 * STDEV  = 1.1  102  ft/sec 

• Projectile  reference  area  (C0DE#  = 122)  is  a 
uniformly  distributed  uncertainty  (TYPE  = 3)  with 
a nominal  or  mean  value  of 

— 3 o 

VALUE  = 3.068  X 10  ftz 
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Eq.  (A. 4-15), formed  from  the 
correlation  coefficients  is 
positive  definite  as 
theoretically  required. 


Figure  3-3  IBM-029  Keypunch  Control  Cards  for  Card  Group 


Table  3-5  Error  Source  Model  Code  Numbers 


ERROR  SOURCE 

STATISTICAL 

DISTRIBUTIONS 

CODE  NUMBER (S) 

•INITIAL  CONDITIONS 
•VELOCITY  VECTOR 

- MAGNITUDE 

- ORIENTATION 
•INERTIAL  ORIENTATION 

- ATTITUDE 

- ATTITUDE  RATE 

GAUSSIAN 

RAYLEIGH/UNIFORM* 

RAYLEIGH/UNIFORM* 

RAYLEIGH/UNIFORM* 

101 

120  & 121 

116  & 117 
118  & 119 

•PHYSICAL  CHARACTERISTICS 

•WEIGHT 

GAUSSIAN 

125 

•MOMENTS  OF  INERTIA 

- AXIAL 

GAUSSIAN 

126 

- PITCH 

GAUSSIAN 

127 

•PHYSICAL  DIMENSIONS 

- REFERENCE  AREA 

UNIFORM 

122 

- LENGTH 

UNIFORM 

123 

- BASE  DIAMETER 

UNIFORM 

124 

•AERODYNAMIC  CHARACTERISTICS 

•STATIC  COEFFICIENTS 

- DRAG  VARIATION 

EFFECTS 

VELOCITY** 

GAUSSIAN 

105  & 107 

ANGLE  OF  ATTACK 

GAUSSIAN 

130 

- NORMAL  FORCE 

GAUSSIAN 

110 

•DYNAMIC  COEFFICIENTS 

- PITCH  DAMPING 

GAUSSIAN 

112 

- MAGNUS  MOMENT 

GAUSSIAN 

113 

•SPIN  RATE 

GAUSSIAN 

128 

•TRIM  ANGLE  OF  ATTACK 

RAYLEIGH/UNIFORM* 

114  & 115 

•STATIC  MARGIN 

GAUSSIAN 

111 

•ATMOSPHERIC  EFFECTS 

•CONSTANT  WINDS 

RAYLEIGH/UNIFORM* 

102  & 103 

* DENSITY  VARIATIONS 

GAUSSIAN 

104 

♦Denotes  a Rayleigh  distribution  of  magnitude  with  a uniform 
360°  distribution  in  orientation. 


♦♦Drag  variation  with  velocity  closely  approximated  by  Cv=  Cv  + 
Kj)/V2  which  is  fit  to  two  data  points  (Cxp  V1)  = (0.03585, 
16.740)  and  (Cx2»  V2)  = (0.11951,  3906.).  Uncertainty  in  Cx^  and 
CX2* 


The  standard  deviation  is 

STDEV  = (2/3%) (3.068  x 10"3)  = 2.045  x 10"  5 ft2 
and 

T0L  ~yjT  * STDEV  = 3.543  x 10“5  ft2 

• Constant  winds  (Code  Numbers  102  and  103)  are 
Rayleigh  distributed  in  magnitude  and  iniformly 
distributed  in  direction.  This  type  of  distri- 
bution is  modeled  as  two  independent  Gaussian 
random  variables  (TYPE  = 2).  There  is  no  average 
wind  so 


VALUE  = 0.0 


As  discussed  in  Appendix  A. 4. 4,  the  equivalent 
Gaussian  standard  deviation  is  given  by 


STDEV  = 


J 

1.2533 


where  J is  the  mean  magnitude  given  in  the  mean 
value  column  of  the  error  source  model,  so 

STDEV  = — ~ — = 8.777  ft/sec 

1.2533 


and 


T0L  - 3 * STDEV  = 2.633  x 10  ft/sec 


These  and  the  other  error  source  Group  2 cards  are  shown  in 
Figure  3-4. 


The  error  source  model  of  Table  1-2  presumed  Fire  Control 
made  no  simplifying  assumptions  in  the  prediction  of  the  pro- 
jectile trajectory.  Thus,  since  Fire  Control  would  always  make 
use  of  the  known  nominal  values,  the  first  Group  2 card  in 
Figure  3-4  sets  a system  control  (C0DE#  = 201)  to  initialize  the 
Fire  Control  trajectory  parameters  at  the  nominal  values  of  the 
Real  World  trajectory  by  declaring  it  a TYPE  = 5 equal  to  1. 

This  considerably  simplifies  the  input.  Even  though  C0DE#  = 201 
denotes  an  integer  variable,  VALUE  is  input  as  a floating  point 
number.  Nominal  range  (C0DE#  = 109)  is  entered  as  a TYPE  = 5 
constant  value,  as  shown  in  Figure  3-4.  The  second  order  drag 
effect  variable  (C0DE#  = 130)  appears  on  two  Group  2 cards. 

The  latter  prevails.  The  dependent  variables  (TYPE  = 7),  force 
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Figure  3-4  Encoded  Error  Source  Model 
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the  statistics  of  dispersion  at  nominal  time  to  be  calculated 
and  printed.  Card  Group  2 is  closed  by  the  required  card  with 
a negative  code  number. 

Card  Group  2 is  now  complete.  The  basic  error  source  model 
has  been  encoded  and  attention  turns  to  selecting  various  HITS 
options.  The  Group  1 card  is  constructed  by  referring  to  Table 
3-1.  An  Analytical  Statistical  mode  calculation  is  requested 
by  the  Group  1 card  in  Figure  3-4.  Card  Group  3 calls  for 
muzzle  velocity  (C0DE#1  = 101)  and  projectile  weight  (CODE#2  = 

125)  variations  to  be  treated  as  negatively  correlated.  The 
Group  3 card  was  constructed  in  accord  with  Table  3-4.  Card 
Group  3 concludes  with  the  mandatory  C0DE#1  = -1  card. 

Exercising  HITS  with  the  case  input  of  Figure  3-4  would 
evaluate  the  total  dispersion  at  nominal  time  due  to  all  twenty 
error  sources  acting  in  consort.  If  it  were  desired  to  evaluate 
the  dispersion  due  to  a single  error  source,  the  case  input 
would  be  generated  by  (1)  duplicating  the  cards  shown  in  Figure 
3-4,  and  (2)  discarding  the  Group  2 cards  pertaining  to  other 
error  sources.  Muzzle  velocity  and  projectile  weight  effects 
would  be  treated  as  a single  source,  since  they  are  declared 
correlated.  An  error  budget  similar  to  Table  1-3  could  be  con- 
structed by  repeating  this  process  for  each  error  source.  The 
resulting  case  inputs  could  be  stacked  and  submitted  simultaneously. 


-60- 


4 . 0 OUTPUT  INTERPRETATION 


This  chapter  presents  four  example  problems  which  demon- 
strate the  basic  modes  of  operation  of  the  HITS  code.  The  user 
is  encouraged  to  reproduce  these  results  to  verify  that  the  code 
has  been  properly  installed  and  to  gain  familiarity. 

4 . 1 Single  Trajectory  Mode 

This  problem  demonstrates  the  evaluation  of  single-shot  dis- 
persion using  the  HITS  code.  This  mode  is  of  interest  to  pro- 
jectile designers  and  exterior  ballisticians . 

Figure  4-1  lists  the  input  cards.  The  first  card  is  the 
Group  1 card  defined  by  Table  3-1.  All  other  cards  belong  to 
Card  Group  2 discussed  in  Table  3-2.  Only  TYPE  = 5 and  TYPE  = 7 
variables  may  appear  in  Single  Trajectory  mode  simulations. 

The  results  of  the  HITS  simulation  are  shown  in  Figure  4-2. 
Frame  'a'  is  the  cover  sheet  which  identifies  the  code  and  sepa- 
rates the  results  of  stacked  cases.  Frames  'b'  through  * f * 
verify  the  inputs.  Frame  'b'  is  the  transcription  of  the  Group  1 
card.  Frame  'c*  lists  the  Card  Group  2 data.  Frames  'd',  'e', 
and  'f'  are  cross  reference  indices  linking  variable  names  and 
code  numbers.  The  objective  is  to  minimize  the  need  to  refer  to 
Table  3-3  for  variable  definitions.  Frame  'd'  gives  the  FORTRAN 
names  and  code  numbers  of  the  dependent  variables  (TYPE  = 7 and 
TYPE  = 8).  Frame  'e'  gives  the  same  information  for  the  deter- 
ministic constants  whose  values  have  been  changed  from  the  presets 
(TYPE  = 5 variables).  All  variables  not  changed  by  the  input 
sequence  are  listed  in  Frame  * f ‘ . 

The  results  of  calculations  resented  in  Figure  4-2g 

and  h.  This  is  the  information  passed  across  the  Trajectory 
Module  interface.  This  level  of  information  may  be  requested  at 
any  time  by  setting  ISPRNT  = 1 in  Card  Group  1 (see  Table  3-1), 
although  it  can  result  in  excessive  output.  The  formats  of 
Frames  'g'  and  'h'  are  identical.  At  the  top  are  the  Fire  Control 
parameters.  These  are  followed  by  the  Real  World  parameters. 
System  Controls,  and  Computed  Quantities.  Table  3-3  gives  the 
definitions  of  all  quantities.  Frame  *g‘  summarizes  the  Fire 
Control  calculation  to  establish  the  coordinates  of  the  nominal 
aim  point  and  the  nominal  time- to- nominal  range.  The  Fire  Con- 
trol summary  is  distinguished  from  the  Real  World  summary  by  the 
value  of  the  IFC  System  Control  (IFC  = 1 — ► Fire  Control, 

IFC  = 0 — ► Real  World).  Frame  'h'  is  the  Real  World  trajectory. 
Values  presented  under  the  Computed  Quantities  heading  quantify 
dispersion. 
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Figure  4-1  Single  Trajectory  Mode  Check  Problem  Input 
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Figure  4-2e  Single  Trajectory  Mode  Check  Problem  Output 
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• Figure  4-2 f Single  Trajectory  Mode  Check  Problem  Output 
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Figure  4-2g  Single  Trajectory  Mode  Check  Problem  Output 
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2h  Single  Trajectory  Mode  Check  Problem  Output 


4 . 2 Range  Check  Mode 


This  problem  demonstrates  the  evaluation  of  dispersion  for 
a systematic  permutation  of  projectile  characteristics.  This 
mode  could  be  used  to  define  a trajectory  as  a function  of  range. 
The  Range  Check  mode  is  of  interest  to  projectile  designers  and 
exterior  ballisticians . 

Figure  4-3  illustrates  the  input  deck  for  a Range  Check 
mode  simulation.  The  first  card  is  the  Group  1 card  defined  by 
Table  3-1.  All  other  cards  belong  to  Card  Group  2 discussed  in 
Table  3-2.  Only  TYPE  = 1,  5,  and  7 variables  may  appear  in  a 
Range  Check  Card  Group  2.  There  must  be  at  least  one  TYPE  = 1 
variable  followed  by  the  Range  Check  values.  There  must  be  at 
least  one  TYPE  = 7 card. 

The  printed  output  generated  by  HITS  is  illustrated  in 
Figure  4-4.  Frame  'a'  is  the  cover  sheet.  Frames  ' b ' through 
* i ' document  the  input.  Frames  'b'  and  'c'  were  discussed  in 
Section  4.1.  Frame  'd'  lists  the  Range  Check  (TYPE  = 1)  vari- 
ables. Frames  'e'  through  'h'  were  discussed  in  Section  4.1. 

As  indicated  in  Frame  'h1,  the  System  Control  IFCRC  is  left  at 
its  preset  value  of  zero.  Thus,  the  Fire  Control  parameters 
remain  at  their  default  values  for  all  Range  Check  combinations 
and  do  not  vary  with  each  combination.  If  a trajectory  were 
being  generated  as  a function  of  range  it  would  be  necessary  to 
set  IFCRC  = 1 (see  Table  3-3  for  details).  Frame  ' i ' lists  the 
Range  Check  combinations.  Since  there  are  two  TYPE  = 1 variables 
with  three  values  each,  HITS  will  evaluate  3*3=9  trajectories. 

Frames  ' j ' and  'k'  summarize  the  results  of  computation. 

The  summary  presents  the  nine  combinations  in  a sequence  of  blocks 
with  the  variable  appearing  last  in  Card  Group  2 varying  most 
rapidly.  The  TYPE  = 1 independent  variables  (IND.  VAR.)  values 
appear  first,  immediately  followed  by  the  TYPE  = 7 dependent 
variable  values.  The  code  numbers  are  linked  to  FORTRAN  names 
by  the  lists  presented  in  Frames  'd'  and  1 e1 . 
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Figure  4-4a  . Range  Check  Mode  Check  Problem  Output 


Figure  4 -4b  Range  Check  Mode  Check  Problem  Output 
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Figure  4-4d  Range  Check  Mode  Check  Problem  Output 
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Figure  4-4g  Range  Check  Mode  Check  Problem  Output 
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Figure  4-4h  Range  Check  Mode  Check  Problem  Output 


Figure  '4-4j  Rarjge  Check  Mode  Check  Problem  Output 


4 . 3 Analytical  Statistical  Mode 


This  problem  demonstrates  the  analytical  evaluation  of 
dispersion  statistics.  This  mode  is  of  interest  to  projectile 
designers  and  those  concerned  with  the  implications  of  projectile 
dispersion  to  overall  weapon  system  effectiveness. 

Pigure  4-5  illustrates  the  input  deck  for  an  Analytical 
Statistical  mode  simulation.  The  first  card  is  the  Group  1 
card  defined  by  Table  3-1.  All  other  cards  belong  to  Group  2 
discussed  in  Table  3-2.  Only  TYPE  = 2,  3,  4,  5,  and  7 cards 
may  appear  in  an  Analytical  Statistical  mode  simulation.  There 
must  be  at  least  one  TYPE  = 2,  3,  or  4 and  one  TYPE  = 7 variable. 
Figure  4-5  illustrates  the  TYPE  = 4 variable  input  format  The 
additional  data  defining  the  probability  density  function  must 
immediately  follow  the  TYPE  = 4 Group  2 card. 

Figure  4-6  illustrates  the  Analytical  Statistical  mode 
printed  output.  Frames  'a'  through  'h'  have  been  described 
previously.  Frames  ' i ' and  * j ' list  the  TYPE  = 4 probability 
density  functions.  Frames  'k'  and  '1'  give  the  Fire  Control 
( IFC  = 1)  and  Real  World  (IFC  = 0)  solutions  for  the  nominal 
conditions.  Frame  'm'  is  a self-explanatory  warning  which 
appears  only  when  code  numbers  320  and  321  are  declared  TYPE  = 7. 

Frames  'n'  and  'o'  of  Figure  4-6  present  the  nominal  case 
independent  and  dependent  variable  values,  respectively.  When- 
ever the  warning  doesn't  appear.  Frame  'o'  is  followed  by  the 
LIMITING  CASES  output  illustrated  in  Figure  4-6,  Frame  'oo', 
presented  at  the  conclusion  of  Figure  4-6.  As  shown  in  Frame 
'oo',  each  independent  variable  is  incrementally  varied  by  the 
T0LERANCE  amount  about  the  nominal  given  by  VALUE.  T0LERANCE 
is  added  to  get  the  first  set  of  ependent  variable  values  and 
subtracted  to  get  the  second.  Return  to  the  example  problem 
at  Frame  'p' . These  are  the  values  of  the  partial  derivatives 
of  the  dependent  variables  with  respect  to  the  independent. 
Independent  variable  code  numbers  are  listed  on  the  left.  De- 
pendent variable  code  numbers  are  listed  above.  1ST  and  2ND 
indicate  the  first  and  second  unmixed  partial  derivatives.  The 
second  order  mixed  partial  derivatives  are  printed  next  as 
illustrated  in  Frame  'q'.  These  are  the  second  order  partial 
derivatives  of  the  dependent  variables  listed  above  with  respect 
to  the  independent  variables  listed  to  the  left.  The  partial 
derivatives  define  the  sensitivity  coefficients  customarily 
listed  in  an  error  budget. 


Frame  * r'  of  Figure  4-6  presents  the  dispersion  statistics. 

This  is  the  culmination  of  the  Analytical  Statistical  mode.  The 
heading  summarizes  the  level  of  the  calculation  as  determined 
by  the  I0PRNT  variable  in  Card  Group  1.  The  mean  values,  variances, 
and  standard  deviations  of  the  dependent  variables  are  printed. 
Standard  deviations  are  printed  a second  time  with  greater  num- 
erical precision. 

Had  any  of  the  independent  variables  been  declared  corre- 
lated by  Card  Group  3,  Frame  ' r'  would  have  been  printed  assuming 
no  correlation.  The  following  two  frames  would  document  the 
effect  of  the  correlation  specified  by  Card  Group  3.  The  dis- 
persion would  be  presented  in  the  format  of  Frame  ' r'.  Printing 
the  dispersion  statistics  with  and  without  correlation  make  it 
possible  to  determine  the  size  of  the  effect  from  a single 
simulation. 


1 4 0 

206  5 0.0000000000 

116  2 0.000000000  3.000000000  1.000000000 

1172  0 .000000000  3.000000000  1.000000000 

125  4 13 

0.004430000  0.01 7520  ODO  0.0S39800O0  0.129500000  0.241900000 

0.3520300000  0.3959000000  0.3520300000  0.2419000000  0.1295000000 

0 . 05398  OODO  0.017520000  0.004430000 

15.0000000- 2  IS. 5000000—2  16.0000000-2  16.5000000-2  17.0000000-2 

17.5000000- 2  18.0000000-2  18.5000000-2  19.0000000-2  19,5000000-2 

20. OOOOOOD— 2 20.5000000-2  21.0000000-2 

123  4 13 

0.004430000  0.01 7S200D0  0.053980000  0.129500000  0.241900000 

0.3520300000  0.3989000000  0.3520300000  0.2419000000  0.1295000000 

0 . 053  98  OCOO  0.017520000  0.004430000 

33.0000000- 2  33.5000000-2  34.0000000-2  34.5000000-2  35.00000D-2 

35.5000000- 2  36.0000000-2  35.5000000-2  3T. 0000000-2  37.5000000-2 

38.0000000- 2  38.5000000-2  39.000000D-2 

105  5 0.035800000 

107  5 0.143400000 

111  5 0,050000000 

320  7 

321  7 

306  7 

307  7 

308  7 

3 09  7 

310  7 

-1 


Figure  4-5  Analytical  Statistical  Mode 
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•Figure  4-6e  Analytical  Statistical  Mode  Check  Problem  Output 


Figure  4-6f  Analytical  Statistical  Mode  Check  Problem  Output 


Figure  4-6g  Analytical  Statistical  Mode  Check  Problem  Output 
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Figure  4-6m  Analytical  Statistical  Mode  Check  Problem  Output 
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Figure  4-6q  Analytical  Statistical  Mode  Chec* 
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Figure  4-6oo  Analytical  Statistical  Mode  Check  Problem  Output 
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This  problem  demonstrates-  the  Monte  Carlo  evaluation  of 
dispersion  statistics  and  probability  distributions.  The  Monte 
Carlo  mode  determines  dispersion  by  performing  numerical  experi- 
ments with  the  full  non-linear  trajectory  equations  rather  than 
the  Taylor  series  approximation  used  by  the  Analytical  Statistical 
mode.  Thus  it  is  inherently  more  accurate,  although  large  num- 
bers of  experiments  would  have  to  be  conducted  to  realize  the 
advantage.  The  number  of  experiments  impacts  cost.  This  mode 
is  of  interest  to  projectile  designers  and  those  interested  in 
projectile  dispersion  statistics. 

Figure  4-7  il3ustrates  the  input  deck  for  a Monte  Carlo  mode 
simulation.  The  first  card  is  the  Group  1 card  defined  by  Table 

3- 1.  The  Group  1 card  sets  MC0PT  = 1,  which  forces  the  Statistical 
Processor  into  the  Monte  carlo  mode  after  using  the  results  of 

the  Analytical  Statistical  mode  to  set  up  the  histogram  boundaries. 
With  the  exception  of  the  first  card,  all  cards  shown  in  Figure 

4- 7  belong  to  Group  2.  Only  TYPE  = 2,  3,  4,  5,  7,  and  8 cards 
may  appear  in  a Monte  Carlo  simulation.  There  must  be  a minimum 
of  one  TYPE  -2,  3,  or  4 and  one  TYPE  = 7 or  8.  Figure  4-7 
illustrates  the  use  of  System  Control  IFCRW  (C0DE#  = 201)  to 
initialize  the  Fire  Control  parameters  at  the  Real  World  nominals. 

Figure  4-8  illustrates  the  printed  output  generated  by  a 
Monte  Carlo  mode  simulation.  Frames  'a'  through  'p'  are  the 
results  of  the  initial  Analytical  Statistical  mode  calculations. 
The  user  is  referred  to  Section  4.3  for  a discussion  of  this 
output.  Subsequent  to  the  Analytical  Statistical  calculations, 
INTERNAL  RANGE  histograms  are  set  up  for  all  TYPE  = 7 and 
TYPE  = 8 variables  with  the  cell  boundaries  determined  by  the 
analytically  computed  mean  values  and  standard  deviations. 

USER  RANGE  histograms  are  also  set  up  for  all  variables  with  the 
cell  boundaries  for  the  TYPE  = 8 variables  determined  by  the 
information  supplied  on  the  Group  2 cards,  as  described  in  Table 
3-2.  At  this  point  the  Monte  Carlo  experiments  begin. 

Figure  4-8,  Frame  'q'  and  Frame  'r'  illustrate  the  summary 
information  printed  for  each  of  the  Monte  Carlo  experiments. 

This  sequence  of  experiments  is  generated  by  using  a random 
number  generator  to  simulate  random  variations  in  the  projectile 
parameters  consistent  with  the  input  statistics.  A given  ran- 
dom sequence  may  be  repeated  by  reusing  the  value  of  IRANN0  on 
the  Group  1 card  (see  Table  3-1  for  instructions).  A different 
saquence  will  result  from  every  selection  of  IRANN0.  The  CASE 
variable  in  Frame  'q'  counts  the  experiments.  The  first  line* 
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labeled  IND.  VAR.,  are  the  stochastic  independent  variable 
values  appearing  in  the  order  given  by  Frame  'd'.  The  line 
labeled  DEP.  VAR.  are  the  corresponding  values  of  the  dependent 
variables  arranged  in  the  order  of  Frame  'e'.  Even  though  only 
a small  number  of  experiments  (100)  were  conducted,  the  summary 
print  required  ten  pages  (only  two  are  included  here).  Thus, 
the  user  should  consider  suppressing  the  summary  print  by  setting 
MCPRNT  = 0 in  Card  Group  1,  whenever  large  numbers  of  experiments 
are  conducted.  The  information  in  the  summary  print  is  used  to 
update  the  histograms  after  each  experiment  but  is  not  stored  in 
any  other  way. 

At  the  conclusion  of  the  Monte  Carlo  experiments,  the  histo- 
gram information  is  printed.  In  general,  there  are  two  rets  of 
histograms.  The  INTERNAL  RANGE  histograms  are  printed  first  and 
are  usually  followed  by  the  USER  RANGE  histograms.  Since  the 
example  problem  input  shown  in  Figure  4-7  contains  no  TYPE  * 8 
variables,  the  sample  output  of  Figure  4-8  contains  only  INTERNAL 
RANGE  histograms.  There  is  no  need  to  document  the  USER  RANGE 
histograms  separately  because  they  are  constructed  in  an  identi- 
cal fashion  and  are  subject  to  the  same  interpretation.  The 
only  difference  is  in  how  the  cell  boundaries  are  set  up,  which 
has  already  been  discussed.  The  histogram  data  for  the  example 
problem  commences  in  Frame  's'  of  Figure  4-8.  The  first  line 
of  print  defines  (1)  the  total  number  of  random  experiments 
performed  (SAMPLE),  (2)  the  number  of  experiments  which  were 
rejected  from  the  INTERNAL  RANGE  histogram  set  because  at  least 
one  of  the  dependent  variable  histogram  ranges  were  exceeded 
(INTERNAL  REJECTS)  and  (3)  the  number  of  experiments  rejected 
from  the  USER  RANGE  histogram  set  for  the  same  reason.  The 
latter  is  zero  in  the  example  because  USER  RANGE  histograms  were 
never  set  up.  The  histograms  making  up  the  INTERNAL  RANGE  his- 
togram set  are  presented  on  a variable-by-variable  basis  in 
Frames  's'  through  'y' . Independent  variable  histograms  pre- 
cede the  dependent  variable  histograms.  Generally,  the  block 
of  data  describing  the  variable,  consists  of  (1)  a line  identi- 
fying the  variable  and  statistics  computed  from  the  histogram, 

(2)  a second  line  defining  the  number  of  experiments  rejected 
from  the  set  because  of  range  limitations  or.  this  variable  and 
Analytical  Statistical  mode  statistics,  and  (3)  multiple  lines 
of  print  defining  the  histogram  on  a cell-by-cell  basis.  The 
columns  define  the  cell  boundaries,  the  number  of  occurrences 
(FREQ),  and  the  probability  density  function  (DIST.  FUNCT.). 
Although  the  format  of  the  block  is  the  same  for  all,  there  are 
some  differences  in  interpretation.  Table  4-1  rigorously  in- 
terprets each  printed  value.  Hash  marks  delimit  the  extent  of 
the  comments  contained  in  the  right  hand  column  of  Table  4-1. 
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Figure  4-7  Monte  Carlo  Mode  Check  Problem  Input 
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Figure  Q-j  Monte  Carlo  Mode  Check  Problem 
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Figure  4-8k  Monte -Carlo  Mode  .Check  Problem  Output 
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Figure  4~8m  Monte  Carlo  Mode  Check  Problem  Output 
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Figure  4-8r  Monte  Carlo  Mode  Check  Problem  Out 
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Figure  4-8s  Monte  Carlo  Mode  Check  Problem  Output 
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^Figure  4-8v  Monte  Carlo  Mode  Check  Problem  Output 
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Figtire  4-8y.  Monte  Carlo  Mode  Check  Problem  Output 


Table  4-1  Monte  Carlo  Mode  Output  Interpretation  (cont'd) 


DESCRIPTION 


5.0 


SUMMARY 


This  report  is  a User's  Manual  for  the  Hypervelocity^  Inf light 
Trajectory  Scatter  (HITS)  computer  code.  The  primary  purpose  of 
the  code  is  to  evaluate  "projectile  dispersion"  which  is  defined 
here  as  the  miss  distance  associated  with  inflight  behavior 
not  anticipated  by  Fire  Control.  Projectile  dispersion  poses  a 
fundamental  limit  on  overall  weapon  system  effectiveness.  The 
code  is  central  to  a scientific,  systematic  approach  to  evalua- 
ting projectile  dispersion.  The  method  consists  of  (1)  construc- 
ting a comprehensive  model  for  each  source  of  projectile  disper- 
sion (i.e.,  the  error  source  model)  and  (2)  processing  it  through 
the  HITS  code  to  obtain  a detailed  projectile  dispersion  error 
budget,  replete  with  sensitivity  coefficients.  The  error  budget 
identifies  the  larger  dispersion  sources  and  those  with  signifi- 
cant potential,  as  well  as  the  total  dispersion.  Both  crossrange 
and  downrange  dispersion  are  computed.  The  error  budget  defines 
dispersion  on  a level  conducive  to  interpretation  and  a compre- 
hensive understanding. 

The  HITS  code  evaluates  projectile  dispersion  statistics 
by  manipulating  the  inputs  to  a set  of  trajectory  equations 
representing  the  true  trajectory  and  comparing  the  computed  flight 
path  to  the  trajectory  predicted  by  a second  set  of  equations 
representing  Fire  Control.  The  true  trajectory  is  varied  in 
accord  with  the  error  source  model;  the  Fire  Control  Trajectory 
is  based  on  nominal  projectile  characteristics.  At  the  selection 
of  the  user,  variations  are  governed  by  either  analytical  or 
Monte  Carlo  statistical  methods. 

The  two  sets  of  trajectory  equations  are  closed  form  approxi- 
mations to  the  full  six  degree  of  freedon  (6  DOF)  equations  of 
motion.  They  consist  of  (1)  a particle  trajectory  with  a velocity 
dependent  drag  coefficient  and  (2)  two  perturbation  equations 
which  correct  for  angle  of  attack  effects.  Extensive  6 DOF  simu- 
lations have  been  conducted  to  verify  the  approximations.  These 
equations  also  have  potential  application  tc  operational  Fire 
Control  computers,  the  computation  of  firing  tables,  and  any  other 
situation  requiring  rapid,  accurate,  low-cost  trajectory 
determination . 


1 "Hypervelocity"  is  used  throughout  this  report  to  indicate  the 
limitation  of  the  present  HITS  code  to  projectiles  whose  velocity 
does  not  become  transonic  at  any  point  along  the  trajectory.  No 
other  connotation  is  intended. 


The  HITS  code  consists  of  two  portions.  One  portion  per- 
forms all  input/output  functions  and  the  statistical  calculations. 
This  portion  is  an  adaptation  of  existing  AVCO  software  (i.e., 
YP-58) . It  brings  to  the  HITS  code  additional  capability  as  a 
projectile  design  tool?  a capability  the  code  would  not  otherwise 
have  had.  The  closed  form  trajectory  equations  constitute  the 
second  portion  of  the  code.  These  equations  were  specifically 
developed  for  dispersion  assessment  under  an  earlier  contract.^- 
The  code  was  assembled  under  the  earlier  contract  and  used  to 
evaluate  the  projectile  dispersion  of  a reference  design.  The 
code,  however,  was  developed  only  to  the  level  of  a research 
program  and  documented  only  to  the  extent  necessary  to  support 
the  analyses  of  the  study.  Detailed  input/output  information 
was  not  presented. 

Under  the  present  effort,  the  HITS  code  was  advanced  to  the 
level  of  a production  code.  The  code  is  maintained  in  the  AVCO 
engineering  computer  code  library  as  Production  Code  5127.  This 
report  is  the  associated  User's  Manual.  It  places  maximum  em- 
phasis on  input/output  information.  The  text  contains  all  the 
day-to-day  necessities  required  to  operate  the  code  and  interpret 
the  results  of  computation.  Input/output  information  is  summar- 
ized by  tables  featuring  a quick-reference  step-by-step  format. 
Example  problems  are  presented  and  discussed.  They  verify  proper 
installation  and  illustrate  the  output  format.  The  appendices 
discuss  theoretical  and  programming  aspects  of  the  code.  They 
contain  the  theoretical  development  of  the  closed  form  trajectory 
equations  and  a complete  listing  of  the  code. 


^Gustafson,  T.  G.,  Crimi,  P.,  Bellaire,  R.  G.,  "Dispersion  of 
Increased  Velocity  Projectiles  - Feasibility  Phase,"  Avco 
Corporation,  AVSD-0200-75-CR,  July  1975. 
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ELEMENTS  OF  PROBABILITY  THEORY  AND  STATISTICS 

This  appendix  presents  an  overview  of  the  pertinent  aspects 
of  probability  theory.  Sections  A.l  through  A. 5 review  basic  con- 
cepts and  theoretical  relationships.  Sections  A. 5 and  A. 6 relate 
the  theory  to  the  practical  problems  associated  with  Monte  Carlo 
methods.  More  detailed  treatments  may  be  found  in  text  books. 1*2 

A.l  Random  Phenomenon  and  Probabilities 


Random  phenomenon  are  those  phenomena  whose  fundamental  pro- 
cesses are  either  not  complex  .ly  understood  or  are  too  complex  to 
define  in  an  entirely  satisfactory  manner.  In  either  case,  the 
scientific  interest  in  the  phenomenon  centers  around  the  "average" 
rather  than  the  "precise"  outcome.  Probability  theory  is  a mathe- 
matical discipline  for  quantifying  random  phenomenon. 

A . 1 . 1 Fundamental  Concepts 

Probability  theory  has  three  building  blocks:  (1)  the  sample 
description  space,  (2)  a collection  of  events,  and  (3)  a probability 
function.  This  section  discusses  these  concepts. 

Sample  Description  Space 

The  sample  description  space  is  the  collection  of  all  possible 
outcomes  of  the  random  phenomena.  Each  occurrence  is  thought  of  as 
the  result  of  an  experiment.  For  instance,  the  sample  description 
space  for  a coin  flipping  experiment  would  consist  of  the  two 
possible  outcomes  (H,  T) . The  sample  description  space  for  a tem- 
perature measurement  taken  at  randomly  selected  time  and  location 
would  include  all  values  from  - co  to  +<»  , {-»  , +®  ). 


1 


1 


I 

1 


Events  I 

■I 

' 1 

Random  events  (or  simply  "Events")  are  the  groupings  of  the  | 

fundamental  outcomes . They  are  subsets  of  the  sample  description  j 

space.  Events  are  defined  by  the  analyst.  For  instance,  a raeteoro-  .! 

logist  interested  in  the  likelihood  of  sub-freezing  temperatures  would 


Iparzen,  E.,  Modern  Probability  Theory  and  its  Applications,  Wiley, 
New  York,  1960. 

2 

Davenport,  W.  B.,  and  Root,  W.  L , An  Introduction  to  the  Theory 
of  Random  Signals  and  Noise,  McGraw-Hill,  New  York,  1958. 


1 


i A 
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select  only  two  events  to  subdivide  the  sample  description  space 
of  all  possible  temperatures:  (-  ® , 32)  and  (32,  too  ). 

Probability  Function 


The  probability  function,  P,  assigns  to  each  event  a number 
between  zero  and  one  corresponding  to  the  relative  likelihood  of 
occurrence  of  the  event.  The  probability  of  an  event,  A,  is  the 

etical  ratio  of  the  number  of  times  it  would  occur,  N,  to  the 
< number  of  experiments,  NT. 


P[A] 


(A. 1-1) 


(Monte  Carlo  procedures  estimate  probabilities  by  counting  the 
results  of  experiments.)  An  event  which  (almost)  never  occurs  has 
a probability  of  zero.  An  event  which  (almost)  always  occurs  has 
a probability  of  one.  For  instance,  if  freezing  temperatures 
occurred  40%  of  the  time  the  probability  of  that  event  would  be 


Pl(— >,  32)]  = 0.40 


(A. 1-2) 


Venn  Diagrams 

A Venn  diagram  is  a graphical  method  of  depicting  the  re- 
lationship between  the  sample  description  space  and  events.  (It 
also  gives  a convenient  method  for  illustrating  the  rules  of 
combination  for  probabilities.)  Figure  A-l  shows  a Venn  diagram. 
There  are  three  events  A,  B,  and  C.  Events  A and  B overlap  so 
A and  B simultaneously  occur  for  some  experimental  outcomes. 

Event  C is  exclusive  of  A and  B,  since  neither  A nor  B can  occur 
when  C occurs . 

A. 1.2  Rules  of  Combination 

Let  the  probabilities  of  the  individual  events  in  Figure 
A-l  be  P(A),  P(B),  and  P(C).  There  are  rules  of  combination. 

Let  A + B denote  the  occurrence  of  either  event,  and  let  AB  de- 
note the  simultaneous  occurrence  of  both  events.  The  probability 
of  A + B is 


P (A  + B)  = P (A)  + P (B)  - P (AB) 


(A. 1-3) 
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Referring  to  Figure  A-l,  A + B is  the  area  encompassed  by  either 
events  A or  B,  and  AB  is  the  area  common  to  both  A and  B.  With 
these  visualizations,  Eq.  (A. 1-3)  simply  states  that  P(A+B)  is 
equal  to  P (A)  + P(B)  less  a correction,  P(AB),  for  the  common 
area  which  is  counted  twice  in  the  sum  P(A)  + P(B). 

Since  probabilities  are  positive  numbers,  Eq.  (A. 1-3)  implies 


p (A  + B)  < P (A)  + P (B) 

(A. 1-4) 

Considering  a cc? lectio  i of  events:  A,  B,  C,  ...  with  the  proba- 
bilities Pf'O,  f(B),  P (c) . . . . , the  probability  that  at  least 
one  of  the  events  will  occur  is  P(A  + B + C + ...)  which  is  re- 
lated to  the  individual  probabilities  through  the  inequality: 


P (A  + B + C + . . • ) <.  P (A)  + P (B)  + P (C)  + . . . (A.  1-5) 


It  is  an  inequality  since  it  may  be  possible  for  at  least  two 
of  the  events  to  occur  concurrently.  If  only  one  of  the  events 
can  occur  the  events  are  said  to  be  mutually  exclusive  and  equality 
holds.  If  at  least  one  of  the  events  must  occur,  the  ensemble  of 
events  is  termed  exhaustive  in  which  case  the  left  side  of  the 
above  inequality  is  identically  one. 

A. 1.3  Conditional  Probabilities 


Considering  two  events,  A and  B,  the  probability  of  their 
simultaneous  or  joint  occurrence  is  P(AB).  The  probability 
P (A/B)  is  the  conditional  probability  of  A given  B and  is  defined 
by: 


P (A/B)  = 


P (AB) 
P (B) 


(A. 1-6) 


provided  P(B)  > 0.  The  conditional  probability  P(A/B),  is  the 
probability  A will  occur  given  that  B has  occurred.  Two  events, 

A and  B,  are  said  to  be  statistically  independent  if  knowledge 
that  B has  occurred  is  of  no  value  in  inferring  that  A will  occur 
so  P(A/B)  = P(A).  This  means 


P(ABC . . . ) = P(A)  P(B)  P(C) . . . 


(A. 1-7) 
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occurrence  is  equal  to  the  product  of  the  individual  probabilities. 
A. 2 Random  Variables 

The  discussion  of  the  previous  section  dealt  wit  a random 
phenomena  in  general.  The  discussion  from  this  point  on  is  re- 
stricted to  random  variables,  i.e.,  those  whose  outcomes  are 
numerically  valued.  The  outside  temperature,  the  instantaneous  ( 
voltage  across  a resistor,  and  the  displacement  of  the  impact  poiht 
from  the  target  point  are  just  three  examples  of  random  variables. _ 

A. 2.1  Distribution  Functions 

A numerically  valued  random  phenomenon  is  one  which  is  com- 
pletely summarized  by  the  value  of  a descriptive  variable  (i.e., 
temperature,  voltage,  or  misdistance  in  the  foregoing  examples). 

The  variable  is  referred  to  as  a random  or  stochastic  variable. 

The  sample  description  space  is  all  real*,  numbers,  and  the  events 


are  all  intervals  and  combinations  of  intervals.  The  probability 
function  is  defined  by  the  (probability)  distribution  function, 
F, 


F (X)  = P(x<X) 


(A. 2-1) 


The  distribution  function  evaluated  at  X is  the  probability  that 
the  random  variable  x has  a value  less  than  or  at  the  most  equal 
to  X.  By  definition,  the  cumulative  distribution  function  has  the 
properties  F(- ® ) = 0,  F(+®  ) = 1 and  F(x)  is  a positive,  non- 
decreasing function  of  x.  (Note:  It  is  customary  to  notationally 
distinguish  between  a random  variable,  x,  and  a value  it  could 
assume,  X,  only  when  demanded  by  clarity.)  A distribution  function 
need  not  be  a continuous  function.  It  may  have  step- type  dis- 
continuities at  discrete  points  corresponding  to  values  with  a 
positive  probability  of  occurrence.  Random  variables  may  be 
"discrete,"  "continuous,"  or  "mixed."  The  remainder  of  this  section 
discusses  these  types. 

A discrete  random  variable  can  take  on  specific  discrete  values 
Xi,  X2»  ...Xn.  The  numerical  value  shown  on  a pair  of  dice  is  an 


example  of  a discrete  random  variable.  Each  experiment  must  result 
in  an  integer  value  between  two  and  twelve.  If  fi  is  the  probability, 
Xi  occurs,  the  distribution  function  is  given  as 


and 


p (x<X) 


(A. 2-2) 


N 

F(XN)  = ^fi  = 1 (A. 2-3) 

i = 1 


If  the  distribution  function  derivative  exists,  the  random 
variable  is  said  to  be  continuous  and  the  (probability)  density 
function  is  defined  according  to 


(A. 2-4) 


The  density  function  at  X,  f(X),  is  the  probability  of  the  event 
that  the  random  variable  x falls  in  the  range 


X < x < X + dx 


(A. 2-5) 


The  probability  of  occurrence  of  any  specific  value  of  a contin- 
uous random  variable  is  zero,  since 


limPlX  < * < X + f) 

f - 0 


f(x)  dx 


(A. 2-6) 


Air  temperature  is  an 
the  likelihood  of  any 
of  the  definitions 


example  of  a continuous  random  variable  since 
specific  temperature  is  zero.  As  a result 


f(x)  dx  = F (oc) 


F (— »°)  = 1 


(A. 2-7) 


The  distribution  of  a mixed  random  variable  contains  both  dis- 
crete values,  with  positive  probability  of  occurrence,  and  a con- 
tinuous distribution. 
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Most  random  variables  treated  by  Monte  Carlo  methods  are 
continuous.  The  density  function  is  estimated  by  dividing  up 
the  range  of  the  random  variable  into  small  intervals  called 
cells.  Histograms  are  then  constructed  which  count  the  number 
of  times  the  value  of  the  random  variable  falls  in  each  cell, 

Ni.  The  exact  value  is  discarded  and  only  the  histogram  is  main- 
tained during  execution  of  the  code.  The  density  function  is 
estimated  as  a piece-wise  constant  function  by  the  relationship 


f(x) 


N; 


Nn 


‘b; 


i+1 


(A. 2-8) 


where  Nt  is  the  total  number  of  trials  and  Xki+i  and  x^  are  the 
upper  and  lower  boundaries  of  the  ith  cell.  It  might  be  argued 
that  Monte  Carlo  methods,  in  effect,  approximate  all  random  vari- 
ables as  discrete  since  only  the  histogram  data  is  saved.  Neverthe- 
less, it  is  most  convenient  to  treat  only  continuous  random  vari- 
ables during  the  theoretical  development.  Equation  (A. 2- 8)  is 
all  that  need  be  recalled  to  make  the  translation  back  to  the 
real  world  of  Monte  Carlo  analysis. 

A. 2. 2 Joint  Distribution  Functions 

It  is  often  necessary  to  describe  two  or  more  random  variables 
in  relation  to  each  other.  These  situations  are  handled  by  treat- 
ing each  of  the  random  variables  xi  as  an  element  of  a random  vec- 
tor x,.  The  joint  (probability)  distribution  function  is  defined 
as 


F(X)  = P (x ! < Xj,  x2  < x2 , «n  £ xN) 


(A. 2-9) 


From  the  definition,  the  following  limits  hold 


lim  F (*)  = 1 

" (A. 2-10) 


n 
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(A. 2-11) 


lim  K(x) 

X | ■ » iv 


*1-1-" 

Xi+i--o 


where  Fi(xi)  is  the  distribution  function  of  the  ith  random  vari- 
able. The  joint  (probability)  density  function  is  defined  as 


a" 

<9x1(9x2...c3xn 


F(x)  = f(s) 


(A. 2-12) 


A. 2. 3 Statistical  Independence 

Random  variables  are  said  to  be  statistically  independent 
whenever  the  joint  distribution  functions  factors  such  that 


or  equivalently 


fu)  - n (*i> 

i = 1 


n 

<«-  n w 

i = 2 


(A. 2-13) 


(A. 2-14) 


where  Fi  is  the  distribution  function  and  fi  is  the  density  function 
of  the  random  variable  xi.  Statistically  independent  random  vari- 
ables are  independent  in  the  sense  that  the  value  of  one  does  not 
affect  the  value  of  another,  and  the  definition  is  equivalent  to  the 
one  given  earlier  for  statistically  independent  random  events. 

A. 3 Averages 

Averages  are  of  great  practical  importance.  These  statistics 
capture  a random  variable's  pertinent  characteristics:  its  average 
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value  and  quantitative  assessments  of  randomness.  Section  A. 3.1 
discusses  the  mathematical  operation  of  expectation,  which  is 
the  means  for  computing  averages.  Specific  averages  (i.e.,  the 
moments)  are  discussed  in  Section  A. 3.2.  Section  A. 3. 3 discusses 
an  important  average  used  to  assess  the  correlation  between  two 
random  variables.  Confidence  limits  are  discussed  in  Section 
A. 3.4. 


A. 3.1  Expectation 


If  g(x)  is  any  function  of  the  random  variable  x,  the  ex-  f 
pectation  (average  value)  of  g(x)  is  denoted  by  e[«(x)]  * For 
a discrete  random  variable 


E[g(x)l  = E 8 (Xj)  f; 

i =>  1 


(A. 3-1) 


where  fi  is  the  probability  of  occurrence  of  X^.  In  the  case  of 
a continuous  random  variable 


E [*<*)] 


J g(x)f(x)dx 

— oe 


(A. 3-2) 


where  f is  the  density  function.  The  development  presented  here 
considers  only  continuous  random  variables.  However,  it  should 
be  kept  in  mind  that  completely  analogous  formulas  exist  for  all 
distribution  types  as  ilJustrated  by  Eq.  (A. 3-1). 

The  expectation  of  a function  of  a random  variable,  as  given 
by  Eq.  (A. 3-2),  is  a weighted  average  of  the  function  values  with 
the  weights  determined  by  the  distribution  of  the  random  variable. 
Expectation  is  a linear  operator,  which  obeys 


E 


TV  hep) 


(A. 3-3) 
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minm-jr.  ■ 


where  ai  is  any  non-random  number, 
multiplication, 


The  distributive  law  of 


II  Et*iM 
1 


(A. 3-4) 


holds  only  if  the  variables  are  statistically  independent. 
The  relation 


E[g(x)]  » g (E  [x]) 


(A. 3-5) 


is  generally  not  true  either. 

The  expectation  of  a vector  valued  function  g(x)  is  defined 
as 


E [g(x)] 


g(x)  f(x)  dx 


(A. 3-6) 


where  f(x)  is  the  joint  probability  density  function  and  integration 
is  with  respect  to  all  variables.  The  expectation  operator  is  a 
linear  operator  with  respect  to  vector  random  variables. 

A . 3 . 2 Moments 

The  quantity  E[*k)  , is  called  the  moment  of  x 

i ao 

xk  £(*)  dx  (A.  3-7) 

In  particular  with  k = 1,  EM'-*  is  the  mean  value  of  the  variable 
x.  The  quantity  defined  by  E[(x-x)k]  is  called  the  ktl1  central 
moment  of  x.  Moments  and  central  moments  are  related.  For  example, 
the  second  central  moment  given  as  e[(x-x)2]  can  expressed  in 
terms  of  the  first  two  momentB 


E[(x-x)2]  = E [x2]  - x2 


(A. 3-8) 


-150- 

■ v...'-.-  . ...-L  ■ 


M aittAiyatMkt 





by  virtue  of  the  fact  that  EM  = * and  the  linear  property  of 
the  expectation  operator. 

Several  of  the  low  order  moments  have  practical  interpre- 
tations. By  definition,  the  first  central  moment  is  zero.  The 
second  central  moment  given  as  E t(x  — x)2 1 is  called  the  variance 
of  x and  is  denoted  a 2 . It  is  a measure  of  dispersion  of  x 
about  the  mean.  The  square  root  of  the  variance  is  called  the 
standard  deviation  and  denoted  by  a 

Based  on  their  frequency  of  general  usage,  the  mean  and 
variance  are  by  far  the  most  important  moments.  Occasionally, 
for  non-zero  mean  random  variables  a coefficient  of  variation  is 
defined 


c 


a 


(A.  3-9) 


Of  less  frequent  usage  are  the  third  and  fourth  central  moments. 

The  third  central  moment  provides  a measure  of  the  asymmetry  of  the 
distribution  about  the  mean  value  and  is  referred  to  as  skewness. 
The  fourth  central  moment  is  referred  to  as  kurtosis  and  provides 
an  additional  measure  of  the  clustering  of  the  distribution  about 
its  mean  value. 

A.  3. 3 Correlation 

Let  x and  y be  random  variables  with  mean  values  x and  y. 

The  expectation  of  the  product 

E[(x  — x)  (y  — y)  ] = 

is  called  the  covariance  of  x and  y,  Cov(x,  y).  If  x and  y are 
statistically  independent,  in  accordance  with  (A. 2-4),  then 


jf  (x  - x)  (y  - y)  f (x,  y)  dx  dy  (A. 3-10) 


Cov(x.y)  = E [x -x]  E [y - yl  = 0 


(A. 3-11) 


The  converse  is  not  true.  That  is,  in  general,  the  vanishing  of 
the  covariance  is  not  sufficient  to  insure  statistical  independence. 
Whenever  the  covariance  is  zero,  the  variables  are  said  to  be 
"uncorrelated. " 
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Frequently  a correlation  coefficient  is  used  rather  than  the 
covariance.  The  correlation  coefficient  is  defined  as 

Cov  (x,  y) 

p*y  = ~~^7y  (A.  3-12) 

where  the  and  oy  are  the  standard  deviations  of  x and  y. 

It  can  be  shown  that  the  correlation  coefficient  always  has  values 
in  the  range  -1  to  1.  Again,  a non-zero  value  of  the  correlation 
coefficient  implies  statistical  dependence  but  a vaninshing  corre- 
lation coefficient  does  not  imply  statical  independence. 

The  interrelationship  between  two  random  phenomena,  x and  y, 
is  describable  to  first  order  by  the  correlation  coefficient,  PXy. 
The  correlation  coefficient  is  most  clearly  understood  by  an  ex- 
ample. One  method  of  assessing  the  presence  of  correlation  between 
phenomenon  y and  phenomenon  x would  be  to  plot  y versus  x.  Any 
clustering  of  data  points  in  the  resulting  "scatter  diagram"  would 
indicate  correlation.  The  correlation  can  be  quantified  by  fit- 
ting a straight  line 

y = ax 

(A. 3-13) 

to  the  data  (assuming  zero  means)  using  least  squares  techniques 
to  determine  a.  The  result  would  be  the  best  linear  prediction 
formula  for  y given  x.  This  formula  is  referred  to  as  the  re- 
gression line.  The  regression  line  formula  can  be  written  in 
terms  of  the  standard  deviations  and  correlation  coefficient 

_y  __  * 

ay=P*y"^  (A.  3-14) 

Equation  (A.3-ls)  shows  that  the  correlation  coefficient  is  the 
slope  of  the  best-fit  straight  line  when  the  data  has  been  nor- 
malized with  respect  to  the  standard  deviations.  Since  both  x and 
y are  random,  the  data  will  be  scattered  about  the  regression  line. 
The  scatter  can  be  quantified  by  calculating  the  root  mean  square 
error 


v i = I 


(A. 3-15) 
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where  i denotes  the  individual  data  points.  The  rms  error  can 
be  shown  to  be 


rms 


(A. 3-16) 


Thus  the  correlation  coefficient  is  a measure  of  the  scatter  about 
the  regression  line.  If  PXy  is  zero,  the  rms  error  is  equal  to 
the  uncertainty  in  the  data  itself,  <Ty.  Thus  uncorrelated  random 
variables  (i.e.,  PXy  = 0)  are  completely  random  relative  to  each 
other  and  no  first  order  interrelationship  exists.  Since  the  rms 
error  must  be  positive,  Eq.  (A. 3-16)  requires 


-1  1 Pxy  * 1 


(A. 3-17) 


If  Pxy  = +1  then  the  rms  error  is  zero,  all  the  data  falls  on 
the  regression  line,  y depends  linearly  on  x,  and  y is  not  random 
relative  to  x. 

A. 3.4  Confidence  Limits 

Confidence  limits  bound  the  uncertainty  in  the  value  of  a 
random  variable.  The  limits  are  selected  to  insure  the  random 
variable  will  fall  inside,  with  a certain  level  of  confidence. 
£ince  a random  variable  x is  expected  to  take  on  its  average  value 
x,  the  K-confidence  limit,  °fc,  is  defined  as  the  value  for  which 


P [ I * - x | <.ak]  = K 


(A. 3-18) 


This  defines  in  terms  of  K and  the  distribution.  The  random 
variable  will  fall  inside  the  interval 


1 • “k  < 1 < * + m 


(A. 3-19) 


with  a relative  frequency  of  K. 
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The  Chebyshev  inequality 


P[|t 


<_  ha]  1 — ■ 


(A. 3-20) 


is  plotted  in  Figure  A-2  along  with  other  representative  distri- 
bution types.  It  can  be  seen  to  be  a conservative  lower  bound. 
Substituting  Eq.  (A. 3-20)  into  (A. 3-18)  and  solving  for  the  con- 
fidence limit  yields 

a 

«k  — . 

V 1 “ K (A.  3-21) 


Thus  a random  variable  always  has  a value  within  2 a of  its  mean  75% 
of  the  time  and  within  3a  of  its  mean  89%  of  the  time.  This 
illustrates  that  the  standard  deviation,  C , quantifies  random- 
ness regardless  of  the  distribution. 

A. 4 Typical  Distributions 

It  is  generally  difficult  to  precisely  determine  the  distri- 
bution of  a random  variable.  However,  it  is  usually  possible  to 
assess  characteristic  properties  which  make  it  possible  to  select 
a distribution  model.  Presented  here  are  four  commonly  used  dis- 
tribution types. 

A, 4.1  Uniformly  Distributed  Random  Variable  (TYPE  =3) 

A uniformly  distributed  random  variable  is  equally  likely  to 
take  on  any  value  in  the  interval  X - T to  3f  + T.  The  density 
function  is 


f(x)  = 


|.-i|<  T 


(A. 4-1) 
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x-x  \±h<r 


The  mean  value  is  x;  T is  referred  to  as  the  tolerance,  and 


a 


T 

V? 


(A. 4-2) 


A uniform  distribution  is  characteristic  of  error  sources 
whose  magnitude  is  limited  by  quality  control  procedures.  Pro- 
jectile length  is  a good  example.  Uncertainties  in  projectile 
length  can  be  easilly  controlled  by  measuring  each  projectile  after 
final  machining.  Those  exceeding  quality  control  limits  would  be 
discarded.  Thus,  off-nominal  lengths  beyond  these  limits  do  not 
occur.  Within  these  limits,  the  length  is  equally  likely  to  be 
any  value. 


A. 4. 2 Gaussian  Distributed  Random  Variable  (TYPE  = 2) 


A Gaussian  distributed  random  variable  is  one  whose  distri 
bution  about  the  mean  has  the  familiar  bell  shape.  The  density 
function  is 


f(x) 


(A. 4-3) 


where  x is  the  mean  value  and  a is  the  standard  deviation.  Figure 
A- 3 illustrates  the  Gaussian  curve.  Since  the  interval  x + 3 <7 
contains  better  than  99  percent  of  the  outcomes,  the  tolerance  T 
is  defined  as 


The  effect  of  a change  in  value  of  a is  indicated  in  Figure  A-4 . 
Increasing  a effects  a simultaneous  lowering  of  the  peak  value  and 
a spreading  of  the  tails  of  the  curve.  Changes  in  the  mean  value 
x results  in  a translation  of  the  curve. 


The  Gaussian  form  is  particularly  interesting  when  sums  of 
random  variables  are  considered.  For  example,  suppose  there  is 
a linear  functional  relation 


z = ax  + by 


(A. 4-5) 


where  a,  b are  constants  and  x,  y are  random  variables.  The 
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FIGURE  A- 3 GAUSSIAN  DENSITY  FUNCTION 


PROBABILITY  DENSITY  FUNCTION 


A 


FIGURE  A-4  A COMPARISON  OF  TWO  GAUSSIAN  DENSITY  FUNCTIONS 
WITH  DIFFERENT  STANDARD  DEVIATIONS 
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variable  z,  by  virtue  of  being  a function  of  random  variables,  is 
also  a random  variable.  If  x and  y are  independent  then  one  can 
conclude  that 

z ■*  a*  + by 

oj  = a2  a 2+  b2*  2 (A. 4-6) 

z x y 


by  application  of  the  principles  previously  presented.  Now  in 
general,  the  knowledge  of  the  mean  value  of  z and  its  variance  ar$ 
insufficient  to  infer  the  distribution  of  z.  For  example.  Figure 
A-5  illustrates  two  different  density  functions  which  have  the 
same  mean  and  variance.  However,  it  can  be  shown  that  if  x and 
y are  Gaussian  distributed  then  z will  also  be  Gaussian.  'Phis 
result  can  be  extended  to  linear  combinations  of  any  number  of 
Gaussian  random  variables. 

The  Central  Limit  Theorem  states  that  under  suitable  conditions 
the  sum  of  arbitrarily  distributed  independent  random  variables 
will  become  Gaussian  as  the  number  of  variables  becomes  large.  The 
necessary  condition  is  that  no  single  term  of  the  sum  can  dominate. 
That  is,  the  variance  of  any  one  term  must  not  be  of  the  same  order 
as  the  sum  of  the  variance  of  the  other  terms.  For  instance,  the 
sum  of  twenty  random  variables,  which  are  uniformly  distributed 
between  0 and  1, 

20 

(A. 4-7) 

i = 1 


is  very  nearly  Gaussian. 

In  summary,  a Gaussian  uncertainty  is  one  whose  distribution 
about  the  mean  has  the  familiar  bell  shape.  It  is  appropriate 
for  effects  formed  from  a large  number  of  independent  random  oc- 
currences. For  instance,  variations  in  muzzle  velocity  are  the 
result  of  a large  number  of  independent  effects  occurring  while 
the  projectile  traverses  the  barrel.  Thus,  the  distribution  of 
muzzle  velocity  would  be  expected  to  be  Gaussian. 

A. 4. 3 Jointly  Gaussian  Random  Variables 

Let  x and  y be  jointly  distributed  Gaussian  random  variables. 
Their  joint  density  function  is 
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PROBABILITY  DENSITY  FUNCTION 


FIGURE  A-5  TWO  PROBABILITY  DENSITY  FUNCTIONS  WITH 
IDENTICAL  MEANS  AND  VARIANCES 


f <x,y) 


l 


(A. 4-8) 


2noxay  (l-pXy)1//2 


Exp 


?2  - 2 px y *y  + y 

2(1  - pxy2) 


cr’2 


where 


(A. 4-9) 


(A. 4-10) 


and  x and  y are  the  mean  values  of  x and  y,  <*x  and  <fy  are  the 
standard  deviations  of  x and  y,  and  PXy  is  the  correlation  co- 
efficient of  x and  y.  Considered  separately,  x and  y are  Gaussian 
distributed.  If  the  correlation  coefficient  is  zero,  the  density 
function  factors  and  x and  y are  statistically  independent.  Thus, 
uncorrelated  Gaussian  random  variables  are  statistically  independent. 

Six  indicies  are  commonly  used  to  summarize  the  randomness  of 
two  jointly  distributed  random  variables.  The  standard  deviation, 
a , has  already  been  described.  There  is  a 50%  probability  the 
point  (x,  y)  will  lie  between  two  parallel  lines  which  are  equi- 
distant from  the  origin  and  are  separated  by  twice  the  Linear  Error 
Probable  (LEP) . The  Circular  Error  Probable  (CEP)  is  the  radius 
of  the  circle  with  a 50%  probability  of  occurrence. 


P [/x2  + y2  <.  C.Ep]  = 0.5 


(A. 4-11) 


The  radium  of  the  80%  circle  is  denoted  R80.  The  mean  radius  J 
is  the  average  displacement  of  the  point  (x,  y)  from  the  origin 


J 


E 


[y?T7 


] 


(A. 4-12) 
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(If  x and  y are  the_cross  range  deflections  expressed  as  a frac- 
tion of  the  range,  J is  referred  to  as  the  "average  jump  angle".) 
The  Radial  Standard  Deviation  (RSD)  is  the  rms  value  of  x and  y 


RSD 


>/ E [x2  + y2  I 


(A. 4-13) 


If  x and  y are  jointly  distributed  Gaussian  random  variables, 
a readily  available  text  book'*-  gives  a detailed  development  of  the 
analytic  expressions  relating  the  CEP  to  <FX,  0y  and  px y.  The  CEP 
is  closely  approximated  by 


CEP  = 0.589  (ox  t-  a ) 


(A.  4-14) 


when  PXy  is  zero  and  <7X  and  Oy  differ  by  no  more  than  80%  of  the 
larger.  Provided  <fx  and  <Fy  = O and  PXy  = 0,  all  four  indices 
are  proportional  to  each  other  as  indicated  in  Figure  A-6.  Each 
of  these  statistics  may  be  used  to  define  a circle  which  have  the 
probabilities  of  occurrence  as  shown  in  Table  A-l. 


The  above  discussion  has  dealt  exclusively  with  two  jointly 
distributed  Gaussian  random  variables.  It  is  possible  to  have  an 
arbitrarily  large  number  of  jointly  distributed  Gaussian  random 
variables:  xj_,  X2,  . . . ,xn.  The  joint  density  function  is  defined 
elsewhere. 2 The  distribution  function  is  completely  determined 
by  the  NxN  covariance  matrix 


Cov  (*) 


o%  2 Cov{xl,x2) 

Cov(x2i*l)  °x22 


(A. 4-15) 


Ipitman,  G.  F.,  Inertial  Guidance,  Wiley,  New  York,  1962 

2Davenport,  N.  B.,  and  Root,  W.  L.,  An  Introduction  to  the  Theory 
of  Random  Signals  and  Noise,  McGraw-Hill,  New  York,  1958. 
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The  covariance  matrix  is  symmetric,  positive  definite  and  com- 
pletely determined  by  N the  standard  deviations  $xi  and  the 
correlation  coefficients  ^XjXj  • 

A. 4. 4 Rayleigh  Distributed  Random  Variable 

Let  J be  a Rayleigh  distributed  random  variable.  The  proba- 
bility density  function  is 


(A. 4-16) 


where  J is  the  mean  value  of  the  random  variable  J.  The  distri- 
bution is  completely  specified  by  J and 


J _L  _ i 7 = 0.52272  j 


(A. 4-17) 


The  Rayleigh  distribution  is  interesting  because  of  its  re- 
lationship to  two  jointly  distributed  Gaussian  random  variables. 
If  x and  y are  identically  distributed  (i.e.,  °x  “ = a ) in- 

dependent (i.e.,  PXy  0)  Gaussian  random  variables,  the  radius 


J - / X2  4 y2 


(A. 4-18) 


ia  Rayleigh  distributed  with  mean  J (i.e,  as  determined  by 
Figure  A-5).  The  orientation  angle 


0 = tan-*  y/* 


(A. 4-19) 


is  uniformly  distributed  between  + r radians.  Thus  the  Rayleigh/ 
uniform  distribution  describes  the  Gaussian  distribution  in  polar 
coordinates . 


A. 5 Random  Variable  Generation 


In  Monte  Carlo  analyses  it  is  necessary  to  generate  random 
variables  with  known  distributions.  It  is  desired  to  make  this 
procedure  as  simple  as  possible  to  conserve  running  time,  simplify 
conversion  to  other  machines,  and  not  overcomplicate  the  programming. 

The  HITS  code  contains  a random  number  generator  which  pro- 
vides uniformly  distributed  random  numbers  in  the  interval  from 
zero  to  one.  This  section  defines  transformations  which  convert 
the  output  of  the  random  number  generator  into  one  of  three  random 
variable  types:  uniformly  distributed,  Gaussian,  or  arbitrary 
( tabular ly  defined)  distributed.  It  is  convenient  to  concurrently 
describe  the  method  if  incrementing  the  histogram  cell  counters 
to  record  the  distribution  of  the  generated  random  sequence.  For 
this  discussion  Nc  is  defined  to  be  the  number  of  histogram  cells. 
Ris  the  uniformly  distributed  number  on  the  interval  zero  to  one 
obtained  from  the  random  number  generator. 

A. 5.1  Uniformly  Distributed  Variables  (TYPE  = 3) 

For  uniformly  distributed  variables  it  is  assumed  x is  the 
mean  value  and  T is  the  tolerance.  The  upper  and  lower  bounds 
defining  the  extent  of  the  histogram  are  given  by 


x0  = 7-  T 


xu  = x + T 


(A. 5-1) 


The  random  variable 


* = 7 + (2R  - 1)T 


(A. 5-2) 


is  uniformly  distributed  on  between  xj  and  xu  . The  histogram 
cell  number  whose  counter  should  be  increased  is 


I -=  Int  [Nc  R]  + 1 


(A. 5-3) 


where  Int  yields  the  largest  integer  which  does  not  exceed  the 
argument. 


A. 5.2  Gaussian  Distributed  Variables  (TYPE  = 2) 


For  the  Gaussian  distributed  variables,  it  is  assumed  x 
is  the  mean  value  and  the  tolerance  T is  3 <r  where  <r  is  the 
standard  deviation.  The  upper  and  lower  bounds  for  the  histogram 
are 


= X - 3 O 
*„"■*  + 3a 


(A. 5-4) 


eo  that  the  total  range  of  the  histograms  is  6 <* 

A uniformly  d„stributed  random  number  R is  converted  into  a 
Gaussian  random  number  by  the  following  procedure  which  uses  the 
conversion  table: 

r^  = 0.0 
i2  - 0.28 
r3  = 0.36 
r4  = 0.41 
r5  « 0.46 

An  index  J is  calculated  according  to 

j = Int  [ 10  R]  r l (A.  5-5) 


r6  = 0.50  rll=  1.00 

r7  = 0.54 
rg  = 0.58 
r9  = 0.64 
r io~  0.78 


where  the  integer  function  Xnt  has  previously  been  defined.  The 
Gaussian  distributed  number  is  then  calculated  from 


rg  : (10  R - J + DOj  + l " rJ>  + rJ 


(A. 5-6) 


Rq  is  Gaussian  with  a mean  of  % and  a standard  deviation  of  1/C. 
The  sample  value,  Rq,  is  then  scaled  to  form  the  sample  of  the 
desired  Gaussian  random  variable 


X = x + (2Rg  - 1)  T 


(A. 5-7) 


and  the  appropriate  histogram  cell  number  is  given  by 


(A. 5-8) 


I = Int[NcRG]  + 1 
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A. 5. 3 Arbitrarily  Distributed  Variable  (TYPE  = 4) 

For  the  arbitrarily  distributed  random  variables , the  proba- 
bility density  function  is  defined  by  a sequence  of  points  xi,  fi. 
It  is  assumed  that  these  values  are  tabulated  with  equal  divisions 
of  the  x coordinate.  Furthermore,  the  average  value  of  the  x co- 
ordinates is  assumed  equal  to  the  mean  value  and  occurs  midway  in 
the  table,  the  probability  of  exceeding  the  x table  is  zero,  and 
the  maximum  value  for  the  distribution  density  function  is  known. 
That  is 

l 

* = Y (X1  + *k)  (A. 5-9) 

1 , , (A. 5-10) 

T - T («t  - «1> 

(fj)max  is  known  (A.  5—11) 


Two  uniformly  distributed  random  numbers  are  generated  and 
are  designated  Ri  and  R2.  Using  Ri,  a tentative  sample  value  is 
generated  from  the  expression 

x = 7 + (2Rt  - l)T  (A. 5-12) 


with  an  associated  histogram  interval  number  being  given  as 


I = Int  [Nc  Rj]  + 1 


(A. 5-13) 


An  index  J is  also  determined  from  the  first  random  number  by  the 
parallel  expression 


J = Int  [ K R j ] + 1 


(A. 5-14) 


where  K is  the  number  of  tabulated  points.  The  index  J is  an 
interpolative  index  which  determines  between  which  two  values, 
xj  and  xj  + 1,  Ri  actually  lies.  Using  this  index  the  value  of 
the  probability  density  function  is  approximated  by  linear  inter- 
polation 


(x-x,) 

f (*)  = h + ^1+1  ~ fy) 

J xj  + 1 ~ xj  } 


(A. 5-15) 


-168 


The  ratio  of  this  value  to  the  maximal  value  is  calculated  and 
compared  with  the  second  random  number.  If  the  relation, 


f 


> 


(A. 5-16) 


is  satisfied  then  the  sample  value  of  the  variable  is  accepted?  if 
it  is  not,  the  sample  value  is  rejected,  two  more  random  numbers 
are  generated  and  the  process  repeated.  This  is  performed  until  an 
acceptable  pair  of  random  numbers  is  found  or  a predetermined  number 
of  trials  (20)  have  been  performed  in  which  case  a system  level 
error  is  generated  and  execution  is  suspended. 

This  method  is  called  Von  Newmann  rejection  sampling.  The 
method  will  generate  samples  which  emulate  the  given  distribution 
function.  However,  it  can  be  inefficient.  Referring  to  Figure  A-7, 
the  area  A^  represents  the  area  under  the  distribution  function  in 
the  transformed  plane  of  the  two  random  numbers.  The  method  consists 
of  accepting  pairs  of  numbers  which  lie  under  the  transformed  dis- 
tribution density  curve. 


FIGURE  A-7  VON  NEWMANN  ^EJECTION  SAMPLING  ILLUSTRATION 


The  area  A2  represents  -the  region  where  rejection  occurs.  A 
sampling  efficiency,  t]  , can  be  defined  as  the  ratio  of  sample 
points  accepted  to  the  total  number  of  trials.  This  efficiency 
is 


7 


Ai 

A i + A2  " Al 


(fiU2T 


(A. 5-17) 
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Typically,  it  can  be  expected  that  on  the  order  of  50%  of  the 
samples  will  be  rejected,  which  is  an  acceptable  rate.  For  very 
peaked  distributions  this  method  could  prove  to  be  inefficient. 

A large  number  of  arbitrarily  distributed  variables  would  com- 
pound the  inefficiency.  If  this  inefficiency  proves  to  be  a 
hinderance,  other  methods  would  be  more  appropriate  and  should 
be  substituted. 

A. 6 Monte  Carlo  Inaccuracies 


This  section  discusses  two  sources  of  error  in  Monte  Carlo 
techniques:  the  grouping  error  and  the  sampling  error.  Guide- 

lines are  developed  for  minimizing  these  errors.  In  summary, 
these  errors  are  small  whenever  the  histogram  cells  are  small  and 
the  sample  size  is  large. 


A. 6.1  Grouping  Errors 

In  r Monte  Carlo  analysis,  the  values  resulting  from  the 
experime  .are  sorted  into  histogram  cells  and  counted  rather 
than  recording  the  precise  values.  The  purpose  is  to  conserve 
computer  storage.  Subsequently,  the  histogram  data  is  used  to 
calculate  moments.  Since  the  exact  values  are  lost  in  the  process, 
an  error,  termed  the  grouping  error,  is  introduced  in  the  moments 


To  illustrate  the  problem,  consider  the  random  variable  x 
with  density  function  f(x).  During  Monte  Carlo  experiments  the 
observed  values  of  x are  grouped  into  cells  which  is  tantamount 
to  assuming  that  all  observed  values  fell  at  the  midpoint  of  the 
cell.  This  is  an  error.  Thus,  in  reality  the  histogram  does  not 
correspond  to  the  true  distribution  of  x but  rather  to  a discrete 
distribution  where  x can  take  only  values  associated  with  the  mid- 
points of  the  cells.  That  is 


X 


i = 1,  . . . , Nc 


(A. 6-1) 


are  the  only  possibilities  where  xj^  are  the  midpoints.  The 
probability  associated  with  each  of  these  values  is 


p 


X < 


f(x)  dx 


(A. 6-2) 


where  xbi  and  xbi+i  are  the  lower  and  upper  cell  boundaries. 
The  ensuing  moment  calculation  based  on  the  histogram  is 
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f(x)  dx 


(A. 6-3) 


X 


n _ 


As  a prelude  to  evaluating  the  grouping  error,  consider  the 
true  nth  moment  of  x given  by  the  integral 


xn  f (x)  dx 


(A.  6-4) 


This  integral  can  be  subdivided  into  the  contributions  from  each 
histogram  cell  to  yield 


f(x)  xn  dx 


(A. 6-5) 


where  x^,  xjj^+i  are  the  boundary  values  of  the  cells  and  it 
is  assumed 


fix)  = o 


< *b_ 


X > 


+ 1 


(A. 6-6) 


The  grouping  error  is  the  difference  between  the  calculated 
moment,  Eq.  (A. 6-3)  , and  the  true  moment  as  given  by  Eq.  (A. 6-5). 
The  exact  magnitude  of  the  grouping  error  in  the  n^  moment  depends 
on  the  specific  distribution.  It  cannot  be  evaluated  in  general 
unless  the  density  function  is  approximated.  A good  first  order 
assumption  is  all  that  is  required:  it  is  assumed  that  the  density 

function,  f(x),  is  constant  over  each  histogram  cell,  i.e., 

-bi  c « < *bi+1  (A. 6-7, 


and  each  of  the  cells  are  assumed  to  be  equal  in  size,  Ax  = xbi+i 
- Xfc£.  Thus,  the  midpoints  of  the  cells  are  given  by 


■a* 


? aasEssasssarea  - a ; 
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(A. 6-8) 


■i  i 


N 

The  grouping  error,  € , is  the  difference  between  the  true  and 
calculated  moments  of  order  N 


N, 


<n  = xn  — xn 


E 

i = i 


r 


(A. 6-9) 


With  the  aid  of  the  identity 


(A. 6-10) 


the  grouping  error  can  be  evaluated  for  low  order  moments.  In 
particular 

c1  = 0 


€ 


2 


(A. 6-11) 


Thus,  to  first  order,  there  is  no  grouping  error  in  the  calculation 
of  the  mean  value;  but,  there  is  an  error  ir  the  calculated  vari- 
ance equal  to  Ax2/12,  where  Ax  is  the  cell  size. 

The  expressions  for  the  grouping  error  given  by  Eq.  (A. 6-11) 
is  quite  accurate.  A more  general  derivation  with  a more  accurate 
approximation  to  the  distribution  function  yields  the  same  result. 
Corrections  for  grouping  errors  are  called  Sheppards  corrections 
and  are  treated  in  greater  detail  in  appropriate  texts  on  statis- 
tics.^- These  corrections  are  not  applied  in  the  present  Monte 
Carlo  code  because  they  can  be  made  arbitrarily  small.  Since 
cell  sizes  are  selected  in  proportion  to  the  standard  deviation 
of  the  variables,  the  percentage  error  due  to  grouping  is  of  the 


leromer,  H. , "Mathematical  Methods  of  Statistics,"  Princeton 
University  Press,  1946. 


order  (10/NC)2  where  Nc  is  the  number  of  histogram  cells.  Thus, 
using  more  than  ten  cells  renders  the  grouping  error  on  the  order 
of  one  percent  of  the  standard  deviation  or  less,  which  is  small 
in  comparison  to  other  sources  of  Monte  Carlo  uncertainties. 

A. 6. 2 Sampling  Errors 

Monte  Carlo  methods  are  exact  only  in  the  limit  as  the  num- 
ber of  experiments  (or  samples)  becomes  arbitrarily  large.  Results 
based  on  finite  sample  sizes  contain  errors  called  sampling  errors, 
because  finite  sequences  of  random  numbers  are  not  representative 
of  the  entire  distribution.  This  section  treats  the  effects  of 
finite  sample  size  on  estimates  of  the  density  function  and  the 
moments.  In  summary,  sampling  errors  can  be  minimized  by  taking 
sufficiently  large  sample  sizes.  Guidelines  are  developed. 

Estimated  Density  Function 

A random  number  generator  seeks  to  generate  random  numbers 
which  are  uniformly  distributed  over  the  interval  zero  to  one. 

Let  a sequence  of  Nt  such  numbers  be  sorted  into  Nc  histogram 
cells  of  equal  size  covering  the  range  zero  to  one.  It  would 
be  expected  that  each  cell  would  occur  Nt/Nc  times,  or  a 
relative  frequency  of  1/NC.  Thus,  for  Nc  = 10  and  N«j  = 100, 

10  values  should  fall  in  each  cell  and  the  indicated  probability 
of  occurrence  of  any  given  cell  would  be  1/10.  However,  when  an 
actual  sequence  of  random  numbers  is  sorted  it  is  generally  found 
that  neither  of  these  theoretical  expectations  is  true.  This 
section  discusses  this  problem.  Only  uniformly  distributed  ran- 
dom variables  are  discussed,  although  the  results  are  generally 
true  regardless  of  distribution  type. 


i 
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The  above  hypothetical  situation  can  be  explained  with  the 
aid  of  the  binomial  distribution.  The  binomial  distribution 
states  that  if  an  event  has  a probability  p of  occurring  and  a 
probability  q = 1 - p of  not  occurring,  then  the  probability  of 
exactly  X occurrences  of  the  event  in  N trials  is 


P(X)  = - 


N! 


X!  (N  - X)! 


pX  qN  — X 


(A.  6-12) 


For  any  particular  cell,  the  probability  that  the  generated  ran- 
dom number  will  fall  in  that  cell  is  p = 1/Nc , with  the  nonoccur- 
rence probability  being  q = 1 - 1/NC.  According  to  Eq.  (A. 6-12), 
the  probability  any  cell  will  occur  exactly  Nt/Nc  times  in  Nt 
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(A. 6-13) 


■where  it  is  assumed  that  N<r/Nc  is  an  integer  quantity.  With  the 
assumption  that  the  number  of  samples  is  large  relative  to  the 
number  of  cells,  Nt  » Nc,  Stirlings  approximation 


N! 


s/lrrli  NN  e-N 


(A. 6-14) 


can  be  used  to  simplify  Eq.  (A. 6-13): 


1 

z 

H 

l 

/ 

LNcJ 

J 2 n Nt  (1  — 1/NC) 

(A. 6-15) 


Using  this  expression  it  is  found  that  for  a sequence  of  100 
random  numbers  sorted  into  10  cells  the  probability  any  cell  oc- 
curs exactly  10  times  is  about  1 chance  in  8.  Increasing  the 
size  of  the  sample  to  1000  doesn't  help,  the  probability  of 
exactly  100  samples  in  a cell  is  about  1 chance  in  24.  Alter- 
nately, sorting  the  sequence  of  100  into  5 cells  doesn't  help 
either;  the  probability  of  any  cell  occurring  exactly  20  times 
is  about  1 chance  in  10.  This  illustrates  that  as  the  sample 
size  is  increased,  or  the  number  of  cells  is  decreased,  the  prob- 
ability of  observing  the  theoretically  expected  number  of  occur- 
rences per  cell  diminishes. 

It  would  appear  there  is  no  means  for  obtaining  a more  exact 
definition  of  the  probability  density  function.  This  is  not  true, 
since  the  density  function  is  more  precisely  defined  in  an  average 
sense  by  larger  sample  sizes,  as  illustrated  by  the  following 
example . 

A sequence  of  100  random  numbers  sorted  into  five  cells  has 
an  expected  number  of  occurrences  of  20  per  cell.  A 5%  deviation 
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I 

in  the  actual  occurrences  per  cell  would  be  either  19,  20  or  21 
I in  a cell.  The  associated  probability  is 


.! 


P(19)  + P(20)  + P (21)  =>  3 P (20)  = 0.3  (A. 6-16) 


Increasing  the  sequence  of  random  numbers  tenfold,  the  expected 
number  of  occurrences  is  200  per  cell  and  a 5%  deviation  would 
span  the  twenty  values  190,  191,...,  210.  The  associated  proba- 
bility is 

210 


20  P (200) 


0.6 


(A. 6-17) 


Thus,  the  chances  of  being  within  5%  of  the  theoretical  are  about 
twice  as  good  with  1000  samples  as  they  are  for  100  samples. 

This  example  illustrates  that  the  number  of  occurrences  tends  to 
the  theoretical  as  the  number  of  trials  tends  to  infinity  in  an 
average  sort  of  way.  Thus,  larger  sample  sizes  improve  the  esti- 
mate of  the  density  function.  A general  rule  of  thumb  is  contained 
in  the  old  saw,  "a  few  hundred  is  too  few  and  a thousand  is  plenty." 

Estimated  Moments 

Monte  Carlo  methods  estimate  the  moments  using  Eq.  (A, 6-3) 
by  approximating  the  density  function  with  the  observed  relative 
frequency 

Ni  l 

f (*)-—• (A. 6-18) 

NT  xbi+1“xbi  1 1 + 1 


so  that 


(A. 6-19) 


where  Ni  is  the  number  of  occurrences  of  the  i^h  cell,  is  the 
total  number  samples,  and  xt>i  and  Xbi+i  are  bhe  lower  and  upper 
bounds  of  the  cell.  Thus,  errors  in  the  estimated  moments  are 
the  compound  effect  of  grouping  errors  and  density  function 
uncertainties . 


Section  A. 6.1  showed  that  grouping  errors  can  be  made  arbi- 
trarily  small  by  using  enough  cells.  This  section  treats  the 
moment  errors  associated  with  uncertainties  in  the  estimated 
density  function.  In  order  to  evaluate  the  errors  in  the  moments 
determined  by  Monte  Carlo  methods,  let 


y = y(*) 


(A. 6-20) 


where  x is  the  independent  random  variable  and  y is  dependent. 
Since  x is  a random  variable  and  y is  a function  of  x,  y is  also 
a random  variable.  Therefore,  y has  a den*ity  function  fy  with 
moments  y,  the  mean  value,  and  • the  variance,  defined  by 


y 


y f(y) 


(A. 6-21) 


(y  - y)2  f(y)  dy 


The  problem  is  to  assess  the  errors  in  the  mean  value  and  variance 
computed  from  a finite  Monte  Carlo  sequence. 


In  applying  the  Monte  Carlo  method  a sequence  of  random 
values:  x^,  X2,  ...,  x^T  is  generated  and  the  resulting  sample 
sequence  y^,  y2,  ...»  yNrp  is  computed.  This  latter  sequence  is 
a sequence  of  random  numbers  which  in  turn  has  a mean 


l 


(A. 6-22) 


If  a second  sample  is  generated  using  a different  sequence  of 
random  numbers,  the  second  sample  will  have  a mean  given  as 


Ny 

F2  = "n7  Yi  (A. 6-23) 


which  in  general  will  be  different  from  the  first.  If  the  process 
is  repeated  many  times  then  a set  of  values,  jx2,  . . . , Mj 


will  result.  This  set  of  values  has  a distribution  function  and 
a mean  value  and  variance  associated  with  it.  The  question 
arises  then  as  whether  the  important  properties  of  the  sample 
mean  distribution  can  be  determined  without  actually  determining 
the  sampling  distribution  itself.  This  would  provide  the  neces- 
sary guidelines  for  minimizing  the  sample  errors  in  the  moments. 
This  can  be  done.  The  sample  mean 


Vi 


(A. 6-24) 


is  a function  of  NT  statistically  independent  variables;  y^,  i ■ 1, 
...,  N.  Forming  the  expectation  of  Eq.  (A. 6-24)  gives  the  mean 
value  of  the  sample  mean 


EUi)  = 


(A. 6-25) 


which,  according  to  Eq.  (A. 6-21),  is 


E[fi]  = y 


(A. 6-26) 


Therefore,  the  mean  value  of  the  sample  mean  is  equal  to  the  mean. 
The  variance  of  the  sample  mean  is  found  similarly.  Thus  noting 


- y = 


(yi  - y) 


squaring  both  sides  and  forming  the  expectation 


Var  (fi)  = E[Qx  -y)2] 


(A. 6-27) 


(A. 6-28) 


Var  (ft) 


E[(yi  - y)2) 


(A. 6-29) 


4n»i.v. 
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Since  E 


(yi  “ y)^  " the  variance  of  the  sample  mean  is 

V 

V»r(fi)  * ~~T. — (A. 6-30) 

Nj 

Thus  for  a sample  size  of  Nrp  of  a random  variable  y,  the  standard 
deviation  of  the  sample  mean  is  ffy/y  where  is  the  standard 
deviation  of  y and  NT  is  the  sample  size.  *0n  the  basis  of  a 
single  Monte  Carlo  sample  of  size  N^,  one  would  report  the  esti- 
mate of  the  mean  value  7 as 


y 


= i* 


(A. 6-31) 


The  °y/  sj Nt  is  termed  a sampling  error  for  the  mean.  To  reduce  the 
sampling  error  by  half  requires  a quadruple  increase  in  the  sample 
size.  By  the  central  limit  theorem,  for  a large  number  of  trials 
the  distribution  of  the  sample  mean  will  tend  to  be  Gaussian. 

It  is  therefore  possible  to  associate  confidence  limits  with  it. 

For  example,  the  true  value  of  y should  lie  within  m + 3 <7y//N^ 
with  95%  confidence.  ' 


Just  as  there  is  a distribution  of  the  sample  mean,  there 
is  also  a distribution  associated  with  the  sample  variance 


s 


2 


(A. 6-32 ) 


The  expected  value  of  the  sample  variance  is 


E[S21  = cry2 


(A. 6-33) 


The  (Nt  - 1)  factor  in  Eq.  (A. 6-32)  is  included  so  the  sample 
variance  is  unbiased;  that  is,  without  an  expected  error.  The 
standard  deviation  of  the  sample  variance  is 


i 


V 


(A. 6-34) 
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i 

To  summarize,  the  moments  computed  from  a Monte  Carlo  sample 
of  size  N<p  have  the  following  characteristics 


y - 


(A. 6-35) 


a 


2 


(A. 6-36) 


where  n and  S are  the  sample  mean  and  sample  variance,  re- 
spectively. The  HITS  code  calculates  m and  for  each  of  the 
variables  of  interest.  The  program  does  not  determine  the  samp- 
ling errors  involved.  The  sampling  errors  are  plotted  in  Figures 
A-8  and  A- 9 for  various  confidence  levels.  Clearly,  samples 
less  than  several  hundred  will  produce  moments  with  substantial 
errors.  Samples  exceeding  a thousand  provide  adequate  accuracy. 

While  the  discussion  of  the  sampling  errors  ir  the  moments 
has  been  based  on  a single  variable,  the  basic  conclusions  are 
applicable  to  the  multiple  variable  case  as  long  as  the  variables 
are  statistically  independent. 
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Section  B.l.x,  an • 

Section  B.1.2. 

B.1.1  Analytical  Statistical  Calculations 

The  theoretical  dispersion  assessment  problem  is  to  statisti- 
cally determine  the  effect  of  independent  variables  on  a dependent 
variable.  For  instance,  the  effect  of  the  projectile  error  source 
model  (independent  variables)  on  the  down-range  dispersion  (de- 
pendent variable).  In  general,  an  exact  theoretical  solution  is 
not  possible.  However,  the  equations  can  be  accurately  approxi- 
mated by  a low  order  Taylor  series  for  reasonable  variations  in  the 
independent  variables,  and  the  mean  value  and  variance  of  the  de- 
pendent variable  can  be  estimated.  The  Analytical  Statistical 
mode  mechanizes  this  concept. 


ERROR  SOURCE 
MODEL 


Taylor  Series  Approximation 


The  relationship  between  the  dependent  variable,  y,  and 
the  independent  variables,  xj_,  X2,  . . „xn,  is  representable  as 
an  algebraic  equation 


y (?)  = y (*i>  x2>  • ■ • > xn) 


(B. 1-1) 


where  x^,  . ,.xn  are  treated  as  the  elements  of  a column  vector 
x.  In  general,  the  equation  is  non-linear.  With  the  assumption 
that  the  function  y is  differentiable  and  reasonably  well  be- 
haved, it  can  be  accurately  approximated  by  a second  order  Taylor 
series1  about  the  mean  value  of  the  independent  variables,  x = 

(xi,  X2 / . < «Xjj)  . 


>y(x) 


y(x) 


dy  _ _ 

+ — — ( x ) (x  - x ) + 
ox  ~ 


+ 


(B.l-2) 


where  T denotes  the  matrix  transpose  and 


dy 

d x 


(x) 


dy 

d x j 


(x)  , 


(B.l-3) 


is  the  gradient  row  vector  of  first  partial  derivatives 


^Sokolnikof f , I.  S.,  and  Sokolnikoff,  E.  S.,  Higher  Mathematics 
for  Engineers  and  Physicists,  McGraw-Hill,  New  York,  1941. 
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and 


(?2y  _ 

— (*) 


J?y 

f)*2  dx  i 


( x ) 


d2 


y 


<9xj  ^*2 


(?) 


(x) 


d2y 

<?Xj  dx 


n 


<x) 


d2y 

dx2  dxn 


( *) 


<?2y 

-5 — r~  (V 

*xn  <>«, 


d2y 

dxn  dx2 


(x) 


(B. 1-4) 


Jh 

d*n2 


(x) 


is  the  Jacobian  Matrix  of  second  partial  derivatives,  both  eval- 
uated at  the  mean  value  of  the  independent  variables,  x. 

Dependent  Variable  Mean  Value 


The  mean  value  is  found  by  talcing  the  expectation  of  Eq. 

(B. 1-2 ) . After  some  manipulation,  the  mean  value  of  the  dependent 
variable  can  be  shown  to  be 


(B.l-5) 


This  equation  states  that  y is  composed  of  the  value  of  y at  the 
mean  of  x plus  second  order  corrections  based  on  the  variances 
and  covariances  of  the  independent  variables.  (Note:  HITS  in- 

cludes the  second  term  whenever  the  control  variable  IOPRNT  ^ 3 
and  the  third  whenever  IOPRNT  ^ 5 . ) 


Dependent  Variable  Variance 


to 


Ate  variance  of  the  dependent  variable  ia  coatputed  according 


oy2  - Ef y2]  - y2 


(B.l-6) 


Equations  (B.l-2)  and  (B. 1-5)  are  substituted  and  manipulated  to 
obtain 


t r ®r 


n — 1 n 


+ 2 


2 S ^<5>  ^ <?)  c°'Ui',i> " 


i » 1 j ■ i+  1 


n — 1 n 


» =.  I t_  * J i « 1 j = i+ 1 L J 


(B. 1-7 ) 


n n - 1 n 


+ 2 


j = 1 i = 1 k = i + 1 L 


Cov  (Xj,  xk)  + 


;£££  [&« 

1 m=i  + l i = 1 k = i+l 


n — 1 ^ n — 1 n 

♦j  L 


d2y  _ 

) _ — («) 

l’*k  d*m 


j = 1 m=) 


Cov  (xj,  xk)  Cov(xj,  xm) 


Equation  (B.l-7)  states  that  is  composed  of  linear  terms  based 

on  the  variance  and  covariances  of  the  independent  variables  and 
a profuse  number  of  second  order  corrections.  (Note:  HITS  in- 

cludes the  first  terra  whenever  IOPRNT  ^ 3 and  the  second  whenever 
IOPRNT  ^ 5.)  The  second  order  corrections  are  exactly  correct 
only  when  the  independent  variables  are  Gaussian,  since  the  de- 
velopment of  Eq.  (B.l-7)  assumed 
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Ejdj  - ijX.j  - i|)(«k  - ik)  (*„-*.) I - 


( ov  (■-,  ij)  im)  i 


♦ Co*  («■,  *1^)  Co»  (*■,  *m)  ♦ Co»  (*j,  *„)  Co*  (*j, 


(B. 1-8) 


(Note:  HITS  includes  the  third  and  fourth  terms  in  Eq.  (B.l-7) 

whenever  IOPRNT  > 6 and  the  fifth  and  sixth  whenever  I OP RUT  > 7.) 

Calculation  of  Derivatives 


In  order  to  estimate  the  mean  and  variance  via  Eqs.  (B.l-5) 
and  (B.l-7),  it  is  necessary  to  evaluate  the  dependent  variable 
and  its  derivatives  at  the  mean  values  of  the  independent  vari- 
ables. The  evaluation  of  y(3c)  requires  one  reference  to  the 
projectile  trajectory  module.  HITS  determines  the  derivatives 
by  manipulating  the  inputs  to  the  projectile  trajectory  module 
in  a systematic  fashion.  A central  difference  schesie  is  used  to 
numerically  approximate  the  derivatives.  The  independent  vari- 
ables are  incremented  one  at  a time  both  positively  and  negatively 
about  their  mean  values.  Letting  Xi  represent  the  independent 
variable  mean  values  and  Ax^  a positive  increment,  two  function 
evaluations,  (i.e.,  calls  to  the  projectile  trajectory  module), 
are  performed  to  give 


y*  5;  ♦ <\*j, . . . , (B.l-9) 

y-  - y ( *,•  *2>  • • • • *i  - A*; *„) 


(Mote:  HITS  sets  Ax^  equal  to  the  input  tolerance  value  T0L  for 

each  independent  variable.)  The  derivatives  are  then  calculated 
according  to 


dy  y+  - y 

d*i  2 /'*xi 


(B. 1-10) 
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(B.l-11) 


d2y 


(D 


f ' - 2y  (?)  + 7 


For  the  second  order  mixed  derivatives,  it  is  necessary  to 
evaluate  the  function  y four  times  for  each  pair  of  independent 
variables.  The  four  function  evaluations  for  the  pair  xi  and  x- 
are 


The  second  mixed  partial  derivative  is 


d2y  _ y++  + y ~ (y+  + y > 

(x)  = — . — — 

dxjdxj  “ Axj  Axj 


(B.l-13) 


Since  the  second  order  mixed  partial  derivatives  are  ^mmetric 
with  respect  to  the  independent  variables,  there  are  n(n-l) 
distinct  derivatives  for  any  dependent  variable,  where  n is  the 
numoer  of  independent  variables. 


ar^sages 
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The  increments  used  in  the  second  order  mixed  partial  deriv- 
ative calculations  are  half  the  values  used  for  the  first  order 
derivatives.  The  reason  is  illustrated  in  Figure  B-2  which 
diagrams  points  in  the  xj_  - Xj  plane.  An  ellipse  can  be  drawn 
through  the  four  points  used  for  the  first  derivative  calculation. 
With  regards  to  the  assumed  Taylor  series  approximation,  the 
ellipse  represents  the  boundary  of  the  region  over  which  the 
linearization  is  valid,  by  virtue  of  the  selection  of  the  Ax's. 
For  consistency,  the  off-axis  coordinates  are  halved  to  insure 
they  will  always  be  interior  to  the  ellipse „ 

Summary 

Equations  (B.l-5)  and  (B.l-7)  are  the  basic  equations  for 
the  Analytical  Statistical  mode  of  the  Statistical  Processor. 

If  there  is  more  than  one  dependent  variable  the  equations  are 
applied  to  each  in  turn.  The  objective  is  to  determine  the  mean 
and  variance  of  all  dependent  variables  from  the  input  statistics 
of  the  independent  variables.  This  forces  the  calculation  of 
first  and  second  order  partial  derivatives,  which  is  a major 
undertaking.  Since  the  Analytical  Statistical  mode  is  exercised 
prior  to  the  Monte  Carlo  mode  (to  set  up  the  histograms),  these 
calculations  are  pertinent  to  the  Monte  Carlo  simulations. 

B.1.2  Monte  Carlo  Calculations 

The  Monte  Carlo  mode  consists  of  three  parts:  (1)  generation 
of  representative  independent  variable  values,  (2)  processing  of 
these  values  to  deter  .line  the  corresponding  dependent  variable 
values,  and  (3)  condensing  the  ensemble  of  solutions  to  a usable 
form.  Representative  sequences  of  the  independent  variables  are 
obtained  from  random  number  generators,  as  discussed  in  Appendix 
A.  The  trajectory  equations  of  Appendix  C relate  the  dependent 
variables  to  the  independent.  This  section  discusses  the  third 
point:  data  condensation  using  histograms. 

Cell  Definition 


In  order  to  achieve  substantia}  condensation,  the  histograms 
must  be  constructed  while  performing  the  Monte  Carlo  experiments. 
In  order  to  do  this,  the  boundaries  of  the  cells  must  be  de- 
termined prior  to  the  experiments.  They  are  determined  'zy  an 
equal  division  of  the  range  of  the  variable  as  defined  by  the 
upper  bound,  xu,  and  the  lower  bound,  x^.  The  i^*1  cell  is 
(xbj_#  xbi+i)«  where  the  boundaries  are  given  by 
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Figure  B-2  Coordinates  for  Derivative  Evaluations 
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where  Nc  is  the  number  of  cells.  For  independent  variables, 
and  xu  are  defined  in  Appendix  A. 5 for  the  various  different 
types.  For  dependent  variables,  (TYPE  = 7),  x^  and  xu  are  de- 
termined by  the  preliminary  Analytical  Statistical  mode  calcu- 
lations. They  are 


1i  = X - 

i 


(B. 1-15 ) 


where  x and  o are  the  calculated  values  of  the  mean  value  and 
standard  deviation.  An  alternate  set  of  values  may  be  input 
(TYPE  =8).  The  cell  number  for  a given  value,  x,  is  determined 
by  evaluating 


(B. 1-16) 


where  Int  is  the  largest  integer  less  than  the  argument.  The 
midpoint  in  each  cell  is  termed  the  "class  mark"  and  iR  denoted 

xmi: 


(B. 1-17 ) 


i 

\ . 
i 

r j 


The  Mo  *te  Carlo  mode  counts  the  number  of  times  each  cell  ;j 
occurs  dv  ing  the  experimental  sequence.  The  probability  density  3 
function  is  estimated  according  to  !j 

■ j 


where  Ni  is  the  number  of  times  cell  “i"  occurred  and 


N 


T 


(B. 1-19) 


is  the  total  of  all  cell  counts  for  this  variable.  Since  it  is 
possible  for  dependent  variables  to  fall  outside  the  range  covered 
by  the  histogram  cells,  (xi,  xu);  is  potentially  less  than  the 
total  number  of  experiments,  Nt. 


Moment  Estimates 


The  sample  mean  and  standard  deviation  are  estimated  from 
the  histograms.  The  mean  value  is 


Nc  » 

f V'  f 1+1 

1 x f(x)  dx  ~ J I xf(x)dx 

-I  „ i r.  1 V 


(B. 1-20) 


which  becomes  the  sample  mean  upon  substitution  of  $q,  (B.l-18), 


N, 


f*  = 


LNi 

Xmi 

i •-=  1 1 


(B. 1-21) 


Using  Eg.  (B.l-17)  the  sample  mean  is 


N„ 


* = Xl  + 


V '~h 


i = 1 


(B.l-22) 


The  variance  is  defined  by 


L 


a2  = | (*-I)2f(x)dx 


(B. 1-23) 
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Upon  substitution  of  Eq.  (B.l-18),  and  using  the  sample  mean 
instead  of  the  mean  value,  Eq.  (B.l-23)  becomes  an  expression 
for  the  sample  variance: 


s2  = 


/ *u  ~ xj 

\ Nc 


tk-n 


i = 1 

B„2  Code  Architecture 


(B. 1-24 ) 


This  section  discusses  the  architecture  of  the  HITS  code. 

The  architecture  was  designed  to  meet  the  requirements  of  the 
Statistical  Processor  computational  procedures  just  described. 

The  fcur  storage  arrays  discussed  in  Section  B.2.1  form  the 
basic  structure.  The  input  processor  of  Figure  B-l  plays  the 
intimate  role  described  in  B.2.2.  Statistical  Processor  functions 
are  presented  in  B.2.3.  Subroutine  definitions  and  linking  are 
detailed  in  B.2.4.  The  facility  for  computer  generated  histo- 
gram plots  is  discussed  in  B.2.5. 

B.2.1  Data  Storage 

The  entire  code  is  designed  around  four  arrays  in  common 
storage.  These  arrays  are 

1)  "OE"  array  - This  one-dimensional  array  is  a common  in- 
put-output storage  area  to  be  shared  by  the  projectile 
trajectory  module  and  the  Statistical  Processor.  That 
is,  all  subroutines  are  written  with  all  their  variables 
equivalenced  to  the  OE  array.  The  variable  "code  numbers" 
are  addresses  in  the  OE  array. 

2)  "IA"  array  - This  bicolumned  array  is  used  to  store 
address  links  between  the  "OE"  array  and  the  additional 
storage  areas  of  the  B and  C arrays  described  next. 

The  first  10  rows  of  the  IA  array  are  reserved  for 
dependent  variables. 

3)  "B"  array  - This  one-dimensional  array  is  used  to  store 
most  input  information  and  some  additional  values. 

4)  "C"  array  - This  one-dimensional  array  is  used  to  store 
some  input  information  and  intermediate  output  infor- 
mation. It  is  a general  scratch  pad  storage  for  the 
Statistical  Processor. 


The  data  transfers  between  the  arrays  and  the  computational 
modules  that  effect  these  transfers  are  indicated  in  Figure  B-3. 

B.2.2  Input  Processor  Functions 

The  Input  Processor  is  a collection  of  subroutines.  It 
determines  the  kind  of  calculations  to  be  performed,  the  sub- 
routines to  be  used,  and  the  data  required.  The  Input  Processor 
reads  the  data,  checks  it  for  completeness  and  takes  appropriate 
action  if  errors  are  detected  or  data  is  found  to  oe  missing. 

If  a complete  set  of  data  has  been  input,  control  is  relinquished 
to  the  Statistical  Processor.  This  section  discusses  the  func- 
tions of  the  Input  Processor  in  detail. 

The  Input  Processor  determines  which  trajectory  module  is 
going  to  be  exercised  (at  present,  there  is  only  the  one  described 
in  Appendix  C,  which  is  specified  by  IY  = '),  clears  the  IA  array, 
and  stores  a list  of  code  numbers  into  the  first  column  of  the 
IA  array  starting  at  row  11.  This  it  a complete  list  of  all  the 
addresses  in  the  OE  array  for  which  data  is  required.  The  first 
10  rows  of  IA  are  reserved  for  variables  in  the  OE  array  which 
are  defined  by  the  input  to  be  dependent  variables.  The  OE 
array  is  initialized  at  preset  (i.e.,  default)  values. 

The  input  cards  are  read  and  processed  one  at  a time. 
Generally,  the  card  contains  a code  number  and  five  values  as 
described  in  the  text.  The  five  values  are  placed  in  the  B 
array  with  the  aid  of  a counter.  The  first  column  of  the  IA 
array  is  searched  to  find  the  row  containing  the  code  number. 

When  it  has  been  located,  the  current  value  of  the  counter  is 
stored  in  the  second  column?  and  the  five  values  are  placed  in 
the  corresponding  location  in  the  B array.  Subsequently,  the 
counter  is  advanced  by  six  and  the  next  data  card  is  read  and 
processed  in  a like  fashion.  The  input  proceeds  until  a card 
having  a negative  code  number  is  encountered  which  signifies 
that  all  data  has  been  read  in.  After  all  the  data  cards  have 
been  read,  each  row  of  the  IA  array  has  two  addresses,  the  first 
is  a location  in  the  OE  array,  (i.e.,  the  code  number),  and  the 
second  is  a location  in  the  B array  where  data  is  stored.  Input 
data  is  extracted  by  searching  the  first  column  of  the  IA  array, 
and  then  using  the  second  column  to  vector  into  the  B array. 

Thus,  the  order  of  the  data  cards  in  immaterial.  Note:  if  any 
rows  of  the  IA  array  are  interchanged  the  information  is  not 
disturbed. 
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If  necessary  data  is  missing,  zeros  will  appear  in  the 
second  column  of  the  IA  array.  The  second  column  of  IA  is 
checked  for  zeros  with  appropriate  action  taken  whenever  a zero 
is  encountered.  At  present,  a preset  value  is  substituted. 

However,  HITS  has  the  facilities  to  access  a "Time  Phased  Data  Base." 
This  option  could  be  activated  if  real-time  interactive  terminals 
are  available. 

This  description  of  the  Input  Processor  is  valid  though 
over  simplified.  Several  additions  and/or  modifications  are 
necessary: 

1)  At  the  time  each  data  card  is  read  and  processed,  the 
variable  value  on  the  card  is  stored  in  both  the  OE 
array  and  B array.  This  saves  a separate  pass  through 
the  data  to  initialize  the  OE  array. 

2)  For  missing  data,  a "Time  Phased  Data  Base"  would  be 
treated  like  a second  source  of  data.  That  is,  when 
the  card  reader  is  exhausted  and  data  is  found  missing, 
the  time  phase  data  base  is  assumed  to  be  implemented 
as  a subroutine  which  can  provide  the  missing  data. 

This  logic  is  present  in  the  program,  but  for  the 
present  all  missing  data  is  filled  in  with  preset 
values . 

3)  Dependent  variable  (TYPE  = 7 or  8)  code  numbers  are 
stored  sequentially  in  the  first  10  rcws  of  the  IA 
array.  For  a TYPE  =7,  nothing  is  stored  in  the  B 
array,  and  the  second  column  of  the  IA  array  is  not 
disturbed.  For  a TYPE  = 8,  the  data  is  stored  as 
previously  described. 

4)  For  a constant  input  value  (TYPE  = 5),  a "-1"  is  in- 
serted in  the  second  column  of  the  proper  row  of  the 
IA  array  to  signify  that  a value  was  read  in.  The 
value  is  then  placed  in  the  OE  array. 

5)  A TYPE  « 1 variable  has  a list  of  values  associated 
with  it.  The  list  is  stored  in  the  C array.  The  TYPE, 
number  of  values,  and  tke  starting  address  for  the  list 
in  the  C array  are  stored  sequentially  in  the  B array. 

The  appropriate  index  of  the  B array  is  stored  in  the 
IA  array,  in  the  appropriate  row  of  the  second  column. 
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Since  no  nominal  value  is  defined  for  this  kind  of 
variable,  the  first  comment  does  not  apply. 

6)  TYPE  8 4 variables  have  a two-dimensional  list  associated 
with  them.  They  are  handled  like  a TYPE  = 1 variable 
(see  comment  five).  However,  after  read  in,  a nominal 
value,  tolerance,  variance,  etc.  are  calculated  and 
these  values  are  stored  in  the  B array  as  though  they 
had  been  read  in. 

7)  Inputs  to  the  projectile  trajectory  module  can  be  either 
floating  point  or  fixed  point.  All  input  values  are 
read  in  floating  point  format.  To  signify  a value  has 
to  be  stored  in  the  OE  array  in  an  integer  format,  the 
negative  of  its  code  number  is  stored  in  the  first  col- 
umn of  the  IA  array.  The  code  number  is  supplied  by 
the  Input  Processor  at  initialization  and  therefore, 
there  is  no  apparent  difference  to  the  user.  Also  note 
that  integers  and  floating  point  numbers  are  stored  in 
the  B and  OE  arrays  with  mixed  formats.  The  method  is 
important  to  understanding  the  details  of  the  code. 

Two  equivalenced  names  are  used,  for  instance,  say 
U (5000)  and  IU(5000,2).  The  first  is  double  precision 
floating  and  the  second  is  single  precision  integer. 

When  storing  an  integer  value  IU  (*,  1)  is  referenced, 
and  when  storing  a floating  number  U (*)  is  referenced. 


B.2.3  Statistical  Processor  Functions 

This  section  describes  the  control  functions  performed  by 
the  Statistical  Processor.  There  are  four  modes  of  operation: 
Single  Trajectory,  Range  Check,  Analytical  Statistical,  and 
Monte  Carlo. 

Single  Trajectory 

After  all  data  has  been  read  by  the  Input  Processor,  the 
nominal  case  is  executed  and  the  results  are  stored  in  the  C array. 
If  the  input  stream  contained  only  TYPE  = 5 (constant)  variables, 
the  nominal  case  values  are  printed  and  a normal  exit  results. 

Range  Che ok 

After  all  data  has  been  read  and  all  missing  values  filled 
in,  the  rows  in  the  IA  array  are  sorted  so  that  TYPE  = 1 (Range 
Check)  variables  occupy  rows  11,  12  and  13.  If  there  are  TYPE  = 1 
variables,  they  are  processed  in  a nested  DO  loop  and  the  values 
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of  the  dependent  variables  for  each  combination  are  stored  in  the 
C array.  After  the  DO  loops  have  been  satisfied,  the  C array 
information  is  printed  and  a normal  exit  results. 

Analytical  Statistical 

If  no  TYPE  = 1 variables  are  input,  the  IA  array  rows  are 
sorted  such  that  all  TYPE  = 2,  3 and  4 variables  are  situated  at 
the  top.  After  this  sort,  the  nominal  case  is  executed.  The 
nominal  case  results  are  stored  in  the  C array. 

Following  the  nominal  case,  limit  cases  are  run  with  the 
independent  variables  sequentially  permuted.  The  resulting  de- 
pendent variable  values  are  stored  in  the  C array  behind  the 
nominal  case  results.  The  derivatives  are  calculated  sequen- 
tially. The  resulting  derivative  values  are  stored  in  the  C 
array  by  overwriting  the  previous  function  evaluations.  Us:lng 
the  derivatives,  the  function  mean  values  and  variances  are  cal- 
culated and  stored  in  the  C array  immediately  after  the  deriva- 
tives. The  derivatives,  means  and  variances  are  then  printed. 

These  processes  can  be  time  consuming  and  require  large 
amounts  of  computer  storage.  For  purposes  of  estimating  the 
storage  and  time  requirements,  let  Np  be  the  number  of  dependent 
variables  (TYPE  * 7 or  8),  and  Ni  be  the  number  of  independent 
variables  (TYPE  = 2,  3 or  4).  As  a general  rule,  storage  re- 
quirements vary  with  the  product  NpNj,  whereas  the  time  con- 
sumption varies  with  Nj.  Table  B-l  states  time  (passes  through 
the  projectile  trajectory  module)  and  storage  requirements.  The 
C array  allocates  a total  of  5,000  locations  to  data  storage  of 
the  outputs  as  well  as  tabular  storage  for  TYPE  * 4 variables. 
This  should  be  more  than  adequate  for  most  purposes.  Should 
more  storage  be  required,  it  is  a straightforward  process  to  in- 
crease the  size  of  the  C array  in  each  subroutine  in  which  it 
appears . 


Table  B-l  Time  and  Storage  Requirements 


Control 

Variable 

Function 

Evaluations 

Required 

Storage  Locations 
(C  Array) 

Nominal  Case 

1 

Np 

I0PRNT  > 1 

Limit  Cases 

2N* 

2NdNj 

I0PPNT  > 2 

Derivatives 

0 

2Nd 

I0PRNT  > 4 

Mixed  Derivatives 

4NI(NI  + 1) 

NdNi(Ni-1)/2 

I0PRNT  > 5 

Covariances 

0 

NI(NI+l)/2 
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Monte  Carlo 


The  Monte  Carlo  mode  performes  four  functions:  It  reads 

additional  input,  verifies  the  input,  conducts  Monte  Carlo  ex- 
periments and  constructs  histograms,  and  outputs  the  results. 

The  input  step  of  the  Monte  Carlo  mode  is  preceded  by  execution 
of  the  Analytical  Statistical  mode.  The  Monte  Carlo  mode  ac- 
quires additional  control  parameters  from  the  input  which  de- 
fine the  sample  size,  number  of  cells,  calculational  option, 
and  allowable  number  of  rejected  trial  solutions. 

The  second  function,  input  verification,  consists  of  re- 
arranging some  of  the  values  generated  by  the  Analytical  Statis- 
tical mode  to  be  compatible  with  the  Monte  Carlo  calculations. 

To  efficiently  describe  the  verification  process,  the  range  of 
all  variables  is  assumed  to  be  known  in  terms  of 

(B. 2-1) 

x — T < x < X + T 


where  x is  the  mean  value  and  T is  the  tolerance.  The  relation- 
ship of  these  variables  to  inputted  and  computed  quantities  has 
been  stated  elsewhere.  These  are  summarized  here  along  with  some 
additions  and/or  modifications.  In  general  for  TYPE  =2,  3 or 
8 variables  the  user  specifies  x,  T,  and  the  variance  cr2,  by 
input.  The  Analytical  Statistical  mode  calculates  these  same 
quantities  for  TYPE  =*  4 and  7 variables.  To  insure  proper 
sampling  and  to  detect  possible  user  errors  the  input  is  veri- 
fied according  to  the  following  rules. 


1)  For  TYPE  = 2 or  7 variables  the  tolerance  T is  taken 
to  be  the  larger  of  the  input  specified  value  or  3 a 
and  the  input  value  of  x is  accepted. 

2)  For  TYPE  = 3 or  8 variables  the  given  tolerance  value 
and  input  value  of  x are  accepted  and  the  variance 
value  is  ignored. 

3)  For  TYPE  = 4 variables  the  midrange  value  and  tolerance 
are  calculated  from  the  extremal  values  of  the  input 
list  with  the  input  values  ignored.  For  example,  if 
ten  values  are  entered 


x - - (x10  + *l) 


1 

T = -J  <x10  “ xl> 


(B.2-2 ) 
(B.2-3) 


-200- 


For  TYPE  = 4 variables  the  maximum  value  of  the  distri- 
bution function  is  also  found  and  stored  in  the  location  re- 
served for  the  calculated  variance  value.  Any  quantities 
changed  in  the  verification  process  are  stored  and  are  subse- 
quently printed  on  the  output.  The  input  verification  step  also 
includes  the  initialization  of  the  histogram  counters  as  well 
as  the  counting  of  the  number  of  TYPE  = 8 variables  to  be 
considered. 

The  third  function  in  the  computational  process  is  the 
conduct  of  the  Monte  Carlo  experiment  and  the  construction  of 
histograms.  This  is  the  heart  of  the  whole  computation.  The 
independent  variable  values  are  generated,  processed  through 
the  projectile  trajectory  module,  and  analyzed  to  increment 
appropriate  histogram  ^ell  counters.  The  process  is  then  re- 
peated over  and  over  until  one  of  the  following  terminating 
conditions  occurs: 

1)  Sampling  difficulties  for  a TYPE  = 4 variable  occur. 

2)  The  number  of  trial  solutions  which  are  outside  of 
the  internally  generated  histogram  ranges  exceeds  a 
user  defined  limit. 

3)  The  number  of  trial  solutions  which  are  outside  of  the 
user  supplied  histogram  ranges  exceeds  a user  defined 
limit. 

4)  The  specified  number  of  experiments  have  been  performed. 

Histograms  are  constructed  for  both  the  dependent  and  inde- 
pendent variaoles  Reasonable  estimates  of  the  ranges  of  all 
variables  are  available  from  input  or  are  internally  computed 
and  cel)  boundaries  are  established.  During  the  execution  of 
the  Monue  Carlo  experiments  it  is  possible  for  one  or  more  de- 
pendent variables  to  fall  outside  the  n nge  covered  by  its 
histogram.  The  results  from  such  an  experiment  should  be  re- 
jected from  all  histograms.  This  requires  the  code  to  take  the 
special  precautions  described  in  the  next  two  paragraphs  when 
incrementing  the  histogram  counters . 

To  accommodate  these  problems  the  following  scheme  is  em- 
ployed. The  basic  histograms  are  stored  in  five  integer  arrays. 
Two  of  these  arrays,  IDl  and  III,  are  used  to  store  the  count 
for  each  cell  belonging  to  the  dependent  and  independent  vari- 
ables, respectively.  For  example,  if  each  variable  range 
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contains  five  cells,  the  counter  for  the  second  cell  of  the 
third  dependent  variable  is  the  twelfth  location  of  IDl.  The 
arrays  XD1  and  111  allow  for  a maximum  of  20  cells  with  a maxi- 
mum of  10  dependent  variables  and  20  independent  variables. 

Two  additional  integer  arrays,  IDT(10)  and  IIT(20),  are  used 
for  temporary  storage  of  appropriate  counter  addresses  in  the 
IDl  and  III  arrays.  As  the  independent  variable  values  are 
generated  for  each  trial  solution,  the  appropriate  counter  ad- 
dresses are  stored  in  the  IIT  array.  After  the  dependent  vari- 
ables have  been  evaluated,  their  counter  addresses  are  stored 
in  the  IDT  array.  If  all  of  dependent  variable  values  are 
within  the  range  of  their  respective  histograms,  the  addresses 
in  the  IDT  and  IIT  arrays  are  used  to  increment  the  appropriate 
counters  in  the  IDl  and  III  arrays.  If  one  or  more  of  the  de- 
pendent variable  values  is  not  within  range,  then  this  trial 
solution  is  discarded  and  none  of  1 3 interval  counters  are 
incremented.  In  this  fashion  the  it-LStograms  count  only  accept- 
able experiments.  A fifth  integer  array,  IRl,  is  used  to  count 
the  frequency  with  which  each  of  the  dependent  variables  causes 
experiments  to  be  discarded. 

The  user  may  desire  histogram  ranges  for  the  dependent 
variables  which  differ  from  those  internally  computed.  To  ac- 
commodate this  possibility  three  additional  integer  arrays, 

ID2,  112  and  IR2,  are  included.  These  arrays  are  replicas  of 
the  basic  arrays,  IDl,  III  and  IRl,  and  are  manipulated  in  ex- 
actly the  same  manner  described  above  with  one  exception.  The 
exception  is  that  iser  specifications  in  the  form  of  TYPE  * 8 
variables  are  use',  instead  of  the  internally  generated  ranges. 
Thus,  if  the  user  specifies  a TYPE  = 8 dependent  variable,  the 
routine  provides  two  sets  of  histograms,  one  set  using  internal 
ranges  and  one  set  using  the  TYPE  * 8 specifications. 

The  fourth  and  final  function  performed  on  the  Monte  Carlo 
mode  consists  of  printing  the  Monte  Carlo  results.  The  histo- 
gram values  are  calculated  and  printed.  The  sample  mean  and 
variance  are  calculated  and  printed  along  with  the  corresponding 
values  from  the  Analytical  Statistical  mode.  The  output  is  on 
a variable  by  variable  basis.  Occasionally  the  program  hangs 
up  in  the  output  section  with  an  indicated  error  of  a divide 
fault.  This  usually  occurs  because  none  of  the  generated  so- 
lutions were  within  the  histogram  ranges.  This  is  usually  due 
to  either  a gross  input  error,  or  an  over  restricted  range  for 
a TYPE  * 8 variable. 


B.2.4  HITS  Subroutines 
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The  HITS  computer  code  consists  of  a MAIN  program  and  sixty- 
two  subroutines.  This  section  briefly  describes  the  subroutines 
and  their  functional  interconnections. 

Table  B-2  liscs  and  briefly  describes  each  of  the  sub- 
routines. Further  details  can  be  gained  from  the  program  list- 
ings of  Appendix  D. 

In  Figure  B-4  is  a diagram  depicting  the  subroutine  linking 
arrangement  of  the  HITS  code.  Referring  to  Figure  b-4,  there 
are  four  decision  blocks  labeled  Dl  through  D4,  each  of  which 
has  several  branch  paths.  The  paths  from  the  decision  blocks 
are  mutually  exclusive  and  whenever  the  code  is  exercised  only 
one  path  is  followed.  Decision  block  Dl  is  used  to  select  the 
projectile  trajectory  module.  At  present,  there  is  only  one 
option  (IY  ® 1);  the  others  are  dummy  facilities.  The  second 
decision  block,  D2,  determines  which  of  the  four  Statistical 
Processor  modes  are  to  be  exercised:  Single  Trajectory,  Range 

Check,  Analytical  Statistical,  or  l <ite  Carlo.  Decision  block 
D2  is  controlled  by  the  TYPE  input  data.  The  decision  block 
D3  is  under  direct  control  of  the  user  and  affects  the  decision 
to  use  a Taylor  series  approximation  of  the  Projectile  Trajectory 
module  for  the  Monte  Carlo  experiments.  Decision  block  D4  de- 
termines whether  separate  histograms  are  to  be  generated  for 
user  specified  histogram  ranges.  This  decision  is  made  on  the 
basis  of  whether  or  not  any  TYPE  - 8 variables  were  input. 

Figure  B-5  diagrams  the  subroutine  linking  arrangement  for 
the  Projectile  Trajectory  module.  Note  the  two  distinct  tra- 
jectory calculations  based  on  the  real  world  and  fire  control 
projectile  characteristics. 

B.2.5  Histogram  Computer  Plots 

The  HITS  computer  code  is  designed  to  facilitate  creation 
of .computer  drawn  histogram  plots.  Since  computer  plot  soft- 
ware is  hardware  specific,  the  present  code  doesn't  do  the  plotting 
Rather  it  gathers  the  data  together  into  one  subroutine  (HYSPLT) 
so  that  autc.natic  plotting  can  be  easily  implemented  by  users  who 
desire  it.  This  section  defines  the  format  of  the  data. 


The  two-dimensional  array  TITLE  (I,J)  contains  the  titles 
for  the  histograms.  The  second  subscript,  J,  denotes  the  variable 


Table  B-2  HITS  Subroutine  List 


SUBROUTINE  FUNCTION 

MAIN  Print  input  data  and  sequence  other  subroutine  calls. 

ARTLU  One-dimensional  table  look  up  routine. 

X 

A2  Computes  J “2  dx  using  Real  World  parameters. 

o 

/* 

a2 dx using  Fire  Control  parameters. 

NAMES  Block  data  routine  to  initialize  output  variable  names. 

CHCKIA  Perform  check  to  see  if  any  data  is  missing. 

CHCKIN  Check  the  input  variables: 

a)  establish  number  of  TYPE  * 8 variables 

b)  force  tolerance  = 3 a for  TYPE  = 2 and  7 
variables 

c)  establish  tolerance  for  TYPE  * 4 variables. 

CNVRT  Convert  uniform  random  numbers  to  Gaussian  random 
numbers  with  the  interval  0 to  1 as  3 o limits. 

CR0SS  Computes  the  Real  World  projectile's  oscillatory 

motion  and  crossrange  (horizontal  and  vertical) 
velocity  and  position  due  to  oscillatory  motion. 

CR0SSF  Computes  the  Fire  Control  projectile’s  oscillatory 
motion  and  crossrange  (horizontal  and  vertical) 
velocity  and  position  due  to  oscillatory  motion. 

CRVFT  Generate  output  variable  values  using  1st  or  2nd 
order  Taylor  series  expansions  of  projectile  tra- 
jectory module. 

DJZfABC  Provide  automatic  Range  Check  calculations . 

D0234  Provide  Analytical  Statistical  mode  calculations. 

EXTRA  Output  routine  for  Range  Check  variable  input  data. 

FC2987  Fire  Control  trajectory  module  interface. 


Table  B-2  HITS  Subroutine  List  (Cont'd) 


SUBROUTINE 


Computes  double  precision  exponents  with  error 
checking. 


FILL* 


FILLIN 


62440 


Calculate  mean  values  and  variances  from  histograms. 

Dummy  routine  to  provide  interface  with  time  phased 
data  base  for  missing  information. 

Computes  fifth  and  sixth  terms  in  Eq.  (B.l-7). 

Computes  third  and  fourth  terms  in  Eq.  (B.l-7). 

Dummy  alternate  trajectory  module  interface. 

Dummy  alternate  trajectory  module  interface. 


G298 


HSTG 


Trajectory  module  interface. 

Increment  histogram  counters  using  internally  estimated 
variable  range  values. 


HSTG2  Increment  histogram  counters  using  user  defined  (TYPE  = 8) 

variable  range  values. 

HYSPLT  Interface  for  automatic  histogram  bar  chart  plots. 


INCARD 


INITIL 


Reads  data  cards  and  decodes  the  TYPE  numbers 

Computes  algebraic  constants  using  Real  World  initial 
conditions . 

Computes  algebraic  constants  using  Fire  Control  initial 
conditions . 

Selects  projectile  trajectory  module,  presets 
data  values  and  initialize  INCARD. 


INLHST  Initialize  all  histogram  and  error  counters  to  zero. 

LEADER  Collects  and  stores  data  for  automatic  histogram  bar 

chart  plotter. 
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SUBROUTINE 


FUNCTION 


Sequence  Monte  Carlo  operations  and  monitor  number 
of  trials  and  number  of  rejects. 

Internally  sorts  input  data  to  move  variables  to  top 
of  storage  arrays. 

Locates  covariances  of  independent  variables  in  the 
C array. 

Locates  mean  values  of  dependent  variables  in  the 
C array. 

Locates  variances  of  dependent  variables  in  the 
C array. 

Locates  variances  of  independent  variables  in  B array. 

Locates  first  order  derivatives  in  C array. 

Locates  second  order  derivatives  in  C array. 

Branch  routine  to  select  projectile  trajectory  module. 

Provides  printed  output  of  the  Monte  Carlo  results. 

Computes  second  term  of  Eq.  (B.l-5)  and  second  term 
of  Eq.  {B. 1-7) . 

Formats  random  numbers . 

Random  number  generator  for  supplying  uniform  random 
numbers  on  interval  0 to  1. 

Real  World  trajectory  module  interface. 

Generates  random  sample- values  for  all  TYPE  = 2,  3 or 
4 variables. 

Formats  and  prints  contents  of  OE  array. 

Prints  standard  deviations  with  greater  precision. 

Data  storage  routine  for  efficient  packing  of  the 
C array. 


SUBROUTINE 


FUNCTION 


TG1795  Dummy  alternate  trajectory  module  preset  routii 

TG2440  Dummy  alternate  trajectory  module  preset  routine. 

TG2987  Initializes  trajectory  parameters  at  preset  (default) 
values . 

TRAJT  Computes  Real  World  projectile  velocity  and  range  as 
a function  of  time,  including  the  influence  of  oscil- 
latory motion  on  drag. 

TRAJX  Computes  Real  World  projectile  velocity  and  time  as 

a function  of  range,  including  the  influence  cf  oscil- 
latory motion  on  drag. 

TRAJXF  Computes  Fire  Control  projectile  velocity  and  time  as 
a function  of  range,  including  the  influence  of  oscil- 
latory motion  on  drag. 

TVX  Computes  projectile  velocity  and  range  as  a function 

of  time  for  a particle  (no  oscillatory  drag)  trajectory. 

T4NTV  Computes  mean  value,  variance  and  tolerance  for 
arbitrarily  distributed  variables. 

VXT  Computes  projectile  velocity  and  time  as  a function 

of  range  for  a particle  (no  oscillatory  drag)  trajectory. 

WIND  Computes  projectile  horizontal  velocity  and  position 

brought  about  by  crosswinds. 

XXC  Computes  the  range  at  which  the  projectile's  oscillatory 

motion  has  converged. 

ZATAN2  Computes  the  arctangent  of  a function,  placing  the 
resulting  angle  in  the  proper  quadrant.  Includes 
provisions  to  avoid  computational  errors  for  the 
special  case  of  non-oscillatory  motion. 

ZZ  Computes  appropriate  value  of  h as  defined  by  Eq.  (C.3-44)] 


associated  with  the  histogram.  These  are,  at  most,  KD  + K234 
variables,  where  KD  is  the  number  of  dependent  variables  (appear- 
ing first)  and  K234  is  the  number  of  independent  variable  (ap- 
pearing last).  The  first  subscript  defines  specific  items  re- 
lating to  the  variable  as  defined  by  Figure  B-6.  The  actual 
histograms  are  stored  in  the  three-dimensional  array  ZDATAD 
(I,J,K).  The  index  K denotes  the  KD  + K234  variables  wi4-  . the 
dependent  variables  appearing  first  and  the  independent  st. 
Figure  B-7  illustrates  the  format  of  a two-dimensional  slice 
of  the  array,  ZDATAD (I, J, *) . The  upper  portion  applies  to  the 
dependent  variables  and  the  lower  defines  the  array  for  inde- 
pendent variables. 


TITLE  (I,*) 


INTERNAL 

NAME 


DESCRIPTION 


Histogram  Indicator 


Code  Number 


X-Axis  Label 


Variable  Type 


Number  of  Rejects 


Re j ects 


Histogram  Mean 


MEAN 


Histogram  Variance 


Estimated  Mean 


(MEAN) 


Estimated  Variance 


(VAR) 


Histogram  Cell  Size 


Interval 


Tolerance 


Number  of  Rejects 


Histogram  Mean 


MEAN 


Histogram  Variance 


Estimated  Mean 


(MEAN) 


Estimated  Variance 


(VAR) 


Histogram  R^nge 


Interval 


Tolerance 


Number  of  Cells 


Maximum  Relative  Frequency 

Indicator  for  Spline 
Smoothing  Option 


1 for  TYPE  - 2, 


aero 
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Figure  B-7 


Histogram  Storage  Array  ZDATAD  (I,J,* 


This  appendix  derives  closed  form  expressions  for  the 
trajectory  of  a hypervelocity  projectile.  The  equations  define 
the  trajectory  in  time  as  well  as  space.  The  closed  form  expres- 
sions result  from  realistic  simplifying  approximations  to  the 
full  six  degree  of  freedom  (6  DOF)  equations  of  motion.  All 
assunq>tions  were  wrified  by  comparison  to  exact  solutions  as 
determined  by  numerical  integration  of  the  equations  of  motion. 

In  all  cases , the  closed  form  equations  agreed  to  within  engi- 
neering tolerances.  The  equations  are  of  general  interest. 

Avco  t\as  employed  them  in:  (1)  the  aerodyr -unic  design  of  a hyper- 

velocir.y  projectile,  (2)  the  interpretation  of  ballistic  range 
test  'iata,  and  (3)  the  evaluation  of  projectile  dispersion. 
Potential  applications  include  incorporation  in  operational  fire 
control  computers,  computation  of  firing  tables,  and  any  other 
situation  requiring  rapid,  accurate,  low-cost  trajectory  deter- 
mination. Section  C.l  presents  background  material  and  an  over- 
view of  the  analytical  development.  Particle  trajectory  equations 
are  presented  in  Section  C.2.  Sections  C.3  and  C.4  present  cross- 
range and  downrange  perturbation  equations,  respectively,  used 
to  correct  the  basic  particle  trajectory. 

J. 1 Introduction 

The  complete  coupled,  nonlinear  equations  of  motion  of 
a projectile  in  free  flight  cannot  be  solved  in  a general  closed 
form.  To  obtain  such  solutions,  it  is  necessary  to  make  simpli- 
fying approximations.  Many  such  solutions  have  been  developed 
*n  the  past  under  various  assumptions  for  application  to  a wide 
variety  of  problems.  These  include  the  gross  evaluation  of 
range-velocity-time  histories,  ballistic  range  data  reduction, 
and  the  development  of  jump  angle  expressions  for  specific  launch 
disturbances.  One  of  the  unique  aspects  of  the  trajectory  equa- 
tions presented  here  is  that  both  downrange  as  well  as  crossrange 
dynamics  are  taken  into  account.  The  key  has  been  to  build  upon 
the  existing  work  and  to  piece  together  various  solutions  in 
such  a manner  as  to  accurately  describe  the  downrange  and  cross- 
range dynamics  of  hypervelocity  projectiles. 


The  trajectory  model  has  three  parts.  The  first  part 
describes  the  basic  particle  trajectory  including  the  effects  of 
constant  velocity  winds  in  both  the  crossrange  and  downrange 


directions.  The  second  part  deals  with  angle  of  attack  oscilla- 
tory motion  and  subsequent  crossrange  deflection  from  the  parti- 
cle trajectory.  The  final  part  of  the  model  makes  use  of  the 
angle  of  attack  history  from  the  oscillatory  motion  solution  to 
compute  range,  velocity,  and  flight  time  perturbations  for  cor- 
rection of  the  basic  particle  trajectory.  These  values  account 
for  the  increase  in  drag  brought  about  by  angle  of  attack  oscil- 
lations. The  three  parts  are  then  superimposed  to  describe  the 
complete  trajectory  history  in  six  degrees  of  freedom. 

The  coordinate  system  used  in  the  trajectory  model  is 
shown  in  Figure  C-l.  A consistent  set  of  symbols  is  used 
throughout  this  appendix  and  is  presented  in  Table  C-l  along  with 
code  numbers  which  relate  the  variables  to  the  HITS  Code. 


C.2  Particle  Trajectory 

This  section  presents  the  basic  particle  trajectory 
closed  form  equations.  The  development  of  the  equations  is  dis- 
cussed in  Section  C.2.1.  Section  C.2. 2 presents  results  which 
verify  the  equations. 


C.2.1  Analytical  Development 


The  particle  trajectory  solution  is  fundamentally  the 
same  as  that  developed  and  reported  in  an  earlier  study,1  with 
the  formulation  extended  to  include  the  effects  of  winds  in  the 
crossrange  and  downrange  directions.  The  primary  assumptions 
involved  in  the  solution  are  the  following: 


• The  drag  coefficient  varies  as  1/V  • 

• The  wind  has  constant  speed  and  direction. 

• The  projectile  weathercocks  into  the  wind  and 
flies  a static  zero  angle  of  attack  trajectory. 

• A uniform  density  atmosphere. 

• The  absence  of  gravity. 


The  drag 
according  to 


coefficient  is  allowed  to  vary  with  velocity 


cD  = cDoo  + 


kd 


(C.2— 1) 


!avco  Systems  Division,  "A  Study  of  the  Feasibility  of  the  Two 
Stage  Light  Gas  Gun  as  an  Air  Defense  Weapon,"  AVSD-0197-73-CR, 
June  1973. 
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SYMBOL 


TABLE  C-l 


LIST  OF  SYMBOLS 


DEFINITION  AMD  COMMENTS 


reference  area,  ft4 


fey 


drag  coefficient 

asymptotic  limit  of  CD  as  V -»  oo 

mean  drag  coefficient  including  oscillatory 
motion  effects 

mean  drag  coefficient  excluding  oscillatory 
motion  effects 


CNa  " ^CD. 


pitching  moment  coefficient  about  center  of 
gravity,  M 

i ' 

-P\1 2 3  AD 

dCy 

pitching  moment  coefficient  slope,  , 


magnus  moment  coefficient. 


, rad" 


, sec 


pitch  damping  derivative. 


1 Denot  es  Neal  World  project  3 Jc  parameter 

2 Denotes  I'ire  Control  projectile  parameter 

3 Denotes  Computed  Quantity 


TABLE  C-l 


LI8T  OF  KVMBOL8  (Cont'd) 


HITS  COOK  NUMBERS 
R.*f.  | F.C.  |C.Q. 


DEFINITION  AND  COMMENTS 


S Cy 

, MC 

ir 


nonul  fore*  coefficient,  1 

7*v2a 


no  10 


normal  force  coefficient  elope,  Cn  , rad”1 

Ifl 


axial  force  coefficient,  "J 

7^v2a 


124  24 


axial  force  coefficient  at  aero  angle  of  attack 

reference  length,  ft 

P**  ,„-l 


206  206 
126  26 


ICNCL 


gravitational  acceleration,  32.174  ft/eec2 
On  v0 

£7  (ICNCL  - 0);  -a^S  (ICNCL  - 1) 

• ^ m 

Aerodynamic  coordinate  frame  selector 

2 

roll  moment  of  inertia,  alug-ft 

2 

pitch  moment  of  inertia,  slug-ft 


-Fv2  AD 
jump  angle,  radians 


TAM1.I:  C-l 

LIST  OF  SYMBOLS  (Cont'd) 


HITS  CODE  NUMDF.KS 


SYMBOL 


DEFINITION  AND  COMMENTS 


130 

30 

471 

472 

123 

23 

411 

416,417 

128 

28 

431 

ill 

11 

441,442, 

443,444 

445 

450 

500 

vertical  component  of  jump  angle,  radians 

horizontul  component  of  jump  angle,  radians 

constant  describing  variation  of  C*  with  angle  of 
attack,  rnd“2 

constant  describing  variation  of  C*  with  velocity, 
aec2/ft2 

parameters  defined  by  Eqs.  (C.3-14)  and  (C.3-15) 

projectile  length,  ft 

71 

, sec 

P V0 

projectile  spin  rate,  rad/sec 

Pcr  projectile  resonant  spin  rate,  rad/sec 

Rl*  r2  parameters  defined  by  Eqs.  (C.3-16)  through 

R3,  R4  (C.3-19) 

Rtrim  static  trim  angle  amplification  factor  due  to  spin 


projectile  static  margin. 


flight  time,  arc 


— , sec 

Vo 


tn  I nominal  flight  time,  sec 


<CP  “ *c* 


perturbation  to  flight  tine  arising  from 
oscillatory  notion,  sac 

velocity,  ft/sec 

downrange  component  of  velocity,  5t/sec 

vertical  component  of  velocity,  ft/sec 

horizontal  component  of  velocity,  ft/sec 

projectile  horizontal  velocity  due  to  croaswind, 
ft/sec 

muzzle  velocity,  ft/sec 

velocity  perturbation  due  to  oscillatory  motion, 
ft/sec 

projectile  weight,  lbs 

downrange  component  of  wind  velocity,  ft/sec 
crossrange  component  of  wind  velocity,  ft/sec 
range,  ft 

range  at  angle  of  attack  convergence,  ft 

projectile  center  of  gravity  location  from  nose,  ft 

projectile  center  of  pressure  location  from  nose, 
ft 

downrange  error,  ft 
vertical  displacement,  ft 


F . C „ 

C.Q  * 

426 

427 

400 

14 

402 

476 

403 

401 

409 

408 

15 

405 

406 

407 

SYMBOL 


DEFINITION  AND  COMMENTS 


component  of  initial  velocity  in  the  Y direction, 
ft/sec 

horizontal  displacement,  ft 

component  of  initial  velocity  in  the  Z direction, 
f t/scc 

projectile  horizontal  displacement  due  to  cross- 
wind,  ft 

pitch  angle  of  attack  in  body  fixed  coordinate 
system,  radians 

static  trim  pitch  angle  of  attack  (zero  spin  rate), 
radians 

trim  pitch  angle  of  attack  when  spinning,  radians 

initial  pitch  angle  of  attack,  radians 

| time  rate  of  change  of  aQ  , rad/sec 

total  angle  of  attack,  radians 

upper  envelope  of  a , radians 

yaw  angle  of  attack  in  body,  fixed  coordinate 
system,  radians 

static  trim  yaw  angle  of  attack  (zero  spin  rate), 
radians 

trim  yaw  angle  of  attack  when  spinning,  radianB 
initial  yaw  angle  of  attack,  radians 
time  rate  of  change  of  p , rad/sec 
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TABLE  C-l 


LIST  OP  SYMBOLS  (Concl'd) 


HITS  CODE  NUMBERS 


SYMBOL 


DEFINITION  AND  COMMENTS 


464,465, 

466,429 

467,463 


437,438, 

439 


460,461, 

462,428 


initial  pitch  orientation,  radians 
initial  pitching  rate,  rad/sec 
A,t  a2,  A a components  of  damping  rates,  sec-3- 

•q. *  *2  phase  angles  defined  by  Eqs.  (C.3-24)  and  (C.3-25), 
respectively,  radians 

p atmospheric  density,  slugs/ft3 

4o>  4i>^2  phase  angles  defined  by  Eqs.  (C.4-15),  (C.3-45), 
and  (C.3-46),  respectively,  radians 

I 

roi  meridianal  orientation  of  static  trim,  radians 

A*  meridianal  shift  of  static  trim  orientation  due  to 
spin,  radians 

*o  initial  yaw  orientation,  radians 

<i>0  initial  yaw  rate,  rac/sec 

*>(,,  oij,  <j2,  A<u  components  of  oscillation  frequency  defined  by 
Eqs.  (C.3-7),  (C. 3-3) , (C.3-4),  and  (C.3-8), 
respectively?  rad/sec 


where  Cd  and  Kjj  are  constants  evaluated  for  the  specific 
projectile  and  velocity  regime.  Knowledge  of  the  projectile 
drag  coefficient  at  two  velocities  within  the  range  of  interest 
allows  determir -tion  of  the  two  constants  Cdqq  and  Kq.  The 
drag  coefficient  approximation  may  also  be  expressed  in  terms 
of  the  ballistic  coefficient,  0 (i.e.,  0 ■ W/CpA),  as  indi- 
cated in  Figure  C-2  which  shows  the  accuracy  of  a typical  fit. 


With  the  inclusion  of  constant  velocity  winds,  the  pre- 
vious^- particle  trajectory  solution  takes  the  form 


V*  - <<Vo 


2fV/'E'CD  (Wxt-*)  Kd  / 21\/BCd  (Wxt- x) 


BCr 


(* 


-)! 


1/2 


+ V. 


(C.2-2) 


t = 


\ tan 


’BCr 


(V0  - wx) 


LV 


v*.-w* 


L vo-"*J 


(C.2-4) 


zw  - <v  - *> 


(C.2-5) 


It  is  most  important  to  notice  that  range  is  the  independent 
variable  in  the  particle  trajectory  equations.  Subsequently, 
velocity  and  time-of-f light  corrections  are  computed  and  applied 
as  a function  of  range  and  not  time.  Equation  (C.2-3)  is  in- 
determinate for  the  case  of  constant  drag  coefficient  (Cj.  ■ CDqq 
and  Kd  » 0).  Numerical  evaluation  of  this  case  can  be  made 
using  any  sufficiently  small  value  for  Kp,  such  that  Kq/V2 < < CD  . 
In  the  presence  of  a wind  component  in  the  downrange  direction, 

Wx,  an  iterative  procedure  must  be  used  to  solve  the  equations, 
since  the  product  of  '.he  unknown  flight  time,  t,  and  Wx  appears 
in  the  expression  for  velocity,  i.e.,  Eq.  (C.2-2). 


^■Avco  Systems  Division,  "A  Study  of  the  Feasibility  of  the  Two 
Stage  Light  Gas  Gun  as  an  Air  Defense  Weapon,"  AVSD-0197-73-CR, 
June  1973. 


C.2.2  Verification 


The  particle  trajectory  equations , including  the  drag 
coefficient  approximation,  were  verified  by  comparing  them  to 
exact  numerical  solutions  of  the  6 DOF  equations  of  motion. 

Figure  C-3  summarizes  one  comparison.  It  shows  the  r ange- velocity- 
time  histories  from  the  analytic  model  and  the  6 DOF  simulation 
agree . 


A separate  comparison  was  made  to  evaluate  the  effects 
of  winds  on  crossrange  error.  This  comparison  (as  well  as  all 
subsequent  check  cases  in  this  appendix)  was  made  using  a smaller 
projectile  whose  aerodynamic  and  mass  properties  are  listed  in 
Table  G-2.  With  a muzzle  velocity  of  11,000  ft/sec  and  a con- 
stant crossrange  wind  velocity,  W2,  of  15  ft/sec,  the  6 DOF  cal- 
culations showed  an  11.2  ft  crossrange  drift  after  2 seconds  of 
flight.  The  analytic  model  agreed  to  within  2 percent. 

C.3  Crossrange  Perturbations 

This  section  presents  the  crossrange  perturbation  equa- 
tions that  are  used  to  correct  the  basic  particle  trajectory. 

The  development  of  the  equations  is  discussed  in  Section  C.3.1. 
Section  C.3.2  presents  results  which  verify  the  equations. 

C.3.1  Analytical  Development 

The  crossrange  perturbation  model  addresses  the  oscilla- 
tory motion  and  corresponding  crossrange  effects.  A solution 
for  the  angle  of  attack  history  was  obtained  from  a NACA  report.1 
This  model  employs  the  following  typical  ballistic  range  assump- 
tions : 

e Constant  velocity 
e Constant  spin  rate 
e Linear  aerodynamics 

e Small  trim  angles  due  to  shape  or  mass  asymmetries 
e Small  angles  (cos  0 * 1,  tan  9 » sin  © *»  0) 
e Uniform  density  atmosphere 
e The  absence  gravity 


^■Nelson, 

Referred 

November 


Robert  L.,  "The  Motion  of  Rolling  Symmetrical  Missiles 
to  a i3ody-Axis  System,"  NACA  Technical  Note  3737, 

1956. 


VELOCITY,  KFT/SEC 


At  first  appearance,  the  constant  velocity  assumption 
may  seem  objectionable  since  significant  velocity  variations 
occur  along  the  actual  trajectory.  The  constant  velocity 
assumption  is  required  to  achieve  a closed  form  solution  with  a 
nonzero  spin  rate.  In  reality,  the  method  of  applying  the  con- 
stant velocity  crossrange  corrections  to  the  particle  trajectory 
minimizes  the  errors  incurred  in  calculating  the  perturbation 
assuming  a constant  velocity.  Velocity  affects  crossrange  dis- 
placement through  its  influence  on  the  aerodynamic  force  in  the 
crossrange  direction.  This  force  is  composed  of  two  components i 
the  drag  force  and  the  normal  force.  For  low  drag,  slender  pro- 
jectiles, the  normal  force  dominates  and  the  drag  force  can  be 
safely  ignored.  The  normal  force  is  proportional  to  the  angle 
of  attack.  Thus,  assuming  small  trim  angles  of  attack,  the 
normal  force  is  only  important  while  the  projectile  is  oscillat- 
ing in  angle  of  attack.  Furthermore,  the  first  few  oscillations 
are  the  most  important  because  subsequent  oscillations  tend  to 
average  out.  Since  these  occur  near  the  muzzle  when  velocity 
is  essentially  constant,  the  constant  velocity  assumption  is 
realistic.  The  crossrange  perturbation  is  applied  to  the  basic 
particle  trajectory  as  a function  of  range.  This  scales  the 
'’cone?  t velocity"  crossrange  perturbation  in  accord  with  the 
downrauge  velocity  variation  of  the  particle  trajectory.  Thus, 
the  crossrange  perturbation  equations  simulate  the  true  velocity 
variation. 

The  solution  given  in  the  NACA  report1  is  in  terms  of 
the  body  fixed  coordinate  system  shown  in  Figure  C-4.  The  form 
of  the  solution  is 


t / A t / 

sin  + vy)  - K2«  ®*n  (<*>2 1 ” — ^ + 


(C. 3-1) 


B m Kje  * cos(<ujt'+  vy)  + K2  e ^ cos  (<02*'  — v?)  + $t 


(C.3-2) 


where  a and  0 are,  respectively,  the  pitch  and  yaw  body 
angles  of  attack  as  illustrated  in  Figure  C-4.  These  equations 
are  written  in  terms  of  the  range  expressed  as  a pseudo  flight 
time,  t1  • X/VQ.  This  effects  the  appropriate  scaling  between 


kelson,  Robert  L.,  "The  Notion  of  Rolling  Symmetrical  Missiles 
Referred  to  a Body-Axis  System,"  NACA  Technical  Note  3737, 
November  1956. 


PROJECTILE  BODY  FIXED  COORDINATE  SYSTEM 


the  crossrange  perturbation  solution  assuming  constant  velocity 
and  the  realistic  velocity  variations  contained  in  the  particle 
trajectory.  The  frequency  and  damping  constants  are  defined  as 
follows : 


where 


W, 


o>  i « <u0  - A<o 


w2  “ wo  + Am 
X J ■ Ag  + AX 


Xj  e Xq  “ AX 


4CMa  / Jx\2  /*««  CMfl\  4CMa  / I,  V 

4 I"  r \ V ) ~ + 1 ' + I “ r \ w 


(C. 3-3) 
(C.3-4) 
(C.3-5) 
(C.3-6) 


2 “i  1/2  \ 1/2 


(C.3-7) 


A co 


A 


o 


AX 


+ 


2a>or 


(C.3-8) 


(C. 3-9) 


(C.3-10) 


where 


CMa  ' -CNa(SM)~  (C.  3-11) 

21y 

I'  - -r  (C.3-12) 

pAD  V2 
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m'rrfriliW 


2 V 


#»A*VC 


(C.3-13) 


The  teras  Kj  and  K2  in  Bqe.  (C.3-1)  end  (c.3-2)  ere  eonstents 
involving  the  initial  conditions 


R2  + R2 


(C.3-14) 
(C. 3-15) 


where 


«k»0t«0  ♦ w2(fl0  - - A2(ao  - atru»rt  ♦ AX  “ **»2(ao  “ "mw*  ~ x2(0o  “ Puim*1 


2{«‘  ♦ A Jr) 


(C.3-16) 


r2 


“ "2^0  ~ “trim5  ~ X2(8o  “ 8tri«^  ♦ AA  ♦ w2(8o  “ ^triJ  ~ X2(oo  ~ “trim*1 


2(w2  ♦ A A2) 


(C.3-17) 


»3  “ " fi*im  ~ R1  (C.  3-18) 

r4  “ ao  ~ atii*  - r2  (C.3-19) 


The  initial  conditions  a Q,  0 0,  a Q,  and  6 0 in  Eqs.  (C.3-16) 
and  (C.3-17)  are  related  to  the  projectile  inertial  orientations 
and  rates  by 


«o  - 


(C.3-20) 


-230- 


r 


0o  - 


v 

o 


(C. 3-21) 


9.  - 


pAg 


'o  ~ “o  2W  Vo  ®o  - p 


(C.3-22) 


Po  * - <Ao  - l^CNaVo0o*P«c 


(C.3-23) 


Hie  phase  angles  and  w2  In  Eqs.  (C.3-1)  and  (C.3-2)  are  com- 
puted according  to 


'•  • *i"'1  (?)  • ““(?) 


(C.3-24) 


i ! 

Ni 


*2 


■ ■i"-‘  (?)  • (?) 


(C.3-25) 


V 

1 


Both  the  sine  and  cosine  definitions  are  given  here  to  establish 
the  proper  quadrants  In  which  the  phase  angles  •'i  and  *2  lie. 

the  terms  “trim  and  pttin  in  Eqs.  (C.3-1),  (C.3-2), 

(C.3-16) , and  (C.3-17)  are  body  fixed  trim  values  for  a spinning 
body  (i.e.,  the  “rolling  trim")  and  are  related  to  the  body 
fixed  static  trim  values,  »st  and  fa  , by 


‘trim  * Rtrim  \/«jT  + 0^.  co*  <*ST  + 


(C. 3-26) 


e 


■ ^fim  + fi2  sin  + A<£) 


triw  "mm 


ST  ST 


(C. 3-27 ) 


where 


^ST  “ 


SID 


-1 


0ST 


/.*  + P2 

ST  PST 


-1 


“ST 


s/a2  + /£ 
ST  ST 


(C.3-28) 
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Fro*  the  work  presented  in  an  AIAA  journal,  ^trin  and  ara 
given  by 


A 0 ■ 


-•[(■- 140’  • ('•4)’]"” 

L/jlYI  l 'Pcr'J 


m ♦ tu 


'•(4) 

['  ■ (4)1 

. r - fe) 

■ • (4) 


[■  ■ (4)'] 


(C.3-30) 


where 


--cMapv|AD-| 

_ 3«y  - y J 


CNo  p a*  CMq  D p A 

4w  8(ly  - Ij) 


~CMc  P A® 

_ 2«y  - I.)  _ 


(C. 3-31) 


(C.3-32) 


^Pettus,  Joseph  J.,  "Persistent  Reentry  Vehicle  Roll  Resonance," 
American  Institute  of  Aeronautics  and  Astronautics  Paper  No. 
66-49,  January  1966. 


The  tarn  R^rim  rapreaenta  the  amplification  factor  fo;.  a apin 
rata,  P,  cloaa  to  tha  critical  (or  raaonant)  frequency,  PCR. 

The  tarn  A a represents  tha  maridianal  shift  of  tha  plana  of  the 
rolling  trim  depending  on  the  value  of  the  spin  rate  with  re- 
spect to  Pcr*  At  resonance  conditions  (P  - PCR) , the  amplifica- 
tion is  maximum  and  the  rolling  trim  plane  is  shifted  by  90° 
with  respect  to  the  plane  of  the  static  trim. 

With  the  angle  of  attack  history  established,  the  lift 
force  is  known.  The  crossrange  displacement  is  found  by  inte- 
gration of  the  equations  of  motion  normal  to  the  flight  path. 

The  remainder  of  this  section  is  devoted  to  this  integration. 

Because  of  the  spinning  motion  of  the  projectile,  the 
body  fixed  coordinate  system  in  which  a and  0 are  described, 
see  Figure  C-4,  rotates  with  respect  to  the  inertial  frame  in 
which  the  trajectory  is  described,  see  Figure  C-l.  At  exit 
from  the  muzzle  (t  “ 0)  the  two  coordinate  systems  are  aligned. 

At  this  instant,  positive  a is  nose  up  and  creates  a lift  force 
in  the  vertical  direction  which  is  the  positive  y direction  in 
the  trajectory  frame.  Owin^  to  the  nature  of  the  body  fixed 
coordinate  system,  positive  0 at  t ■ 0 is  nose  left  and  produces 
a horizontal  yaw  force  in  the  negative  z direction.  With  posi- 
tive spin  being  "right  wing  down, " the  body  fixed  lift  and  yaw 
forces  rotate  with  respect  to  the  trajectory  coordinate  frame 
as  illustrated  in  Figure  C-5.  The  force  components  in  che  iner- 
tial trajectory  coordinate  system  are  written  as  follows. 3- 


CN  PA  V' 

(1 

i . a 

a 2 


(C. 3-33) 


CNa  PA 


(C.3-34) 


4>  - p t' 


(C. 3-35) 


Fy  - La  co»^  + sini£ 


(C. 3-36) 


Fx  - La  *in<£  - Lpco*<)> 


(C. 3-37) 


1 In  most  applications,  it  is  more  appropriate  to  use  C^a  instead 
of  Cfla  in  the  ensuing  development.  This  affects  only  the 
definition  of  h,  as  noted  in  Eq.  (C.3-44) . 


Thus  the  equations  of  motion  are 


w d2Y 
* dt'2 


CNa  PA  V2 

[a  co«  P t ' + 0 sin  P t '] 


(C.3-38) 


i»  d2z 

8 dt'2 


CNfl  P A V2 

[a  sinPt'  - ^cosPt'] 


(C.3-39) 


where  a and  0 are  given  by  Eqs.  (C.3-1)  and  (C.3-2).  The  nota- 
tion t'  is  maintained  in  Eqs.  (C.3-38)  and  (C.3-39)  since  it  is 
not  the  true  flight  time. 


Integrating  Eqs.  (C.3-38)  and  (C.3-39)  gives  the  solu- 
tions for  crossrange  velocity  and  displacement  brought  about  by 
oscillatory  motion 
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C 31 
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h dt,l{ 


(C. 3-40) 
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(C.3-41) 


| xj  + («1  + p): 


Ait  r 

* 1 sin  [(a>i  + P)  t'  + vi  + 2^j] 


A*  + (a>2  - Pr 


e ^ tin  [(ciij  — P)t'  — uj  + 2^1 


Kj  sin  (i/j  + 0j)  *2  *in  (^2  — ^ 1 dy 
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where 
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(C. 3-46) 


TI  KIT  numerical  verification  results  presented  in  this  appendix 
use  Cjj  . This  is  a worst  case  assumption  for  the  reference 
projectile. 
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again  take  note  of  the  respective  quadrants  for  and  <t>  2 . 

The  "jump  angle"  is  defined  as  the  root  sum  of  the 
squares  of  the  coefficients  of  the  linear  term  in  Eqs.  (C.3-42) 
and  (C.3-43)  divided  by  the  velocity.  It  is 

JA  - V) Ay  + Ja2  (C.3-47) 

where 


K2  »in  (^2  ■“  ^2) 

>/a  2 + (<^2  — 


Kj  sin  (i/j  + 

+ (ct)  j + P)^ 


(C. 3-48) 
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J 4-  (cjj  4-  P)^ 


4- 


K2  cos  (^2  “ v2> 

\j\\  + (w2  - P)2 


(C. 3-49) 


C.3.2  Verification 


Verification  of  the  crossrange  perturbation  equations 
was  established  by  comparison  to  four  6 DOF  simulations.  The 
results  are  summarized  in  Table  C-3.  In  all  cases,  the  analytic 
solution  agrees  to  within  6%.  An  additional  comparison  is  shown 
in  Figure  C-6,  in  which  crossrange  error  due  to  trim  angle  of 
attack  is  shown  as  a function  of  spin  rate.  The  analytic  model 
is  seen  to  be  in  excellent  agreement  with  the  6 DOF  calculations. 

C.4  Downrange  Perturbations 

This  section  presents  the  downrange  perturbation  equa- 
tions that  are  used  to  correct  the  basic  particle  trajectory. 

The  development  Of  the  equations  is  discussed  in  Section  C.4.1. 
Section  C.4. 2 discusses  results  which  verify  the  equations. 

C.4.1  Analytical  Development 

The  downrange  perturbation  equations  model  the  down- 
range  effects  of  angle  of  attack  oscillations.  In  addition  to 
the  assumptions  stated  in  Sectior  C.3.1,  the  following  assump- 
tions are  required  to  achieve  an  accurate  closed  form  solution: 


13,760  FT 


LIFT 


• Trim  angles  of  attack  cause  negligible  drag. 

• Downrange  effects  can  be  captured  by  mean 
drag  coefficients. 

• Pitch  oscillations  are  lightly  damped. 

In  addition  to  \;.hc  crossrange  perturbation,  oscillatory 
motion  has  a downrange  effect.  Angle  of  attack  oscillations 
cause  the  projectile  -co  experience  a net  increase  in  drag.  The 
exact  mechanism  is  illustrated  in  Figure  C-7.  As  shown  in 
Figure  C-7 (a),  drag  is  defined  as  the  force  acting  in  line  with 
the  velocity  vector.  There  are  two  contributions  when  the  body 
is  at  angle  of  attack,  a component  of  the  normal  force,  CN*  acts 
in  the  drag  direction,  increasing  the  drag  force.  The  second 
component  is  the  axial  force  coefficient,  C*,  which  in  yeneral 
is  a function  of  angle  of  attack,  as  illustrated  in  Figure  C-7(b). 
The  drag  coefficient  is 


Cq  = C%  cos  a + sin  a 


(C.4-1) 


which,  for  small  angles,  becomes 


CD  ^ C*  + CN a 


(C.4-2) 


Modeling  the  axial  force  as  parabolic  in  angle  of  attack, 
cx  - cx  + kc  a2  , and  the  normal  force  as  linear  in  angle  of 
attack  (i.e.,  cN  = cN  a ) the  drag  coefficient  becomes 

Cq  = cx  + (CN  + Kf~  )a 2 (C.4-3) 

o 'a  x 


The  term  is  the  value  of  the  drag  coefficient  at  zero  angle 
of  attack.  It  is  given  by  the  drag  model  of  Section  C.2.1, 

Eg.  (C.2-1).  For  the  general  case  in  which  oscillations  occur 
in  both  the  pitch  and  yaw  planes,  the  angle  of  attack  in  these 
expressions  becomes  the  total  angle  of  attack  (i.e.,  the  angle 
between  the  body  axis  and  the  velocity  vector),  denoted  here- 
after as  cr  . For  small  angles 


a2  = a2  + ft2 


(C.4-4) 
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FIGURE C -7  DOWNRANGE  FORCES  DUE 
TO  ANGLE  OF  ATTACK 


where  a and  0 are,  respectively,  the  pitch  and  yaw  angles  of 
attack  in  the  body  fixed  coordinate  system.  Upon  substitution 
of  Egs.  (C.3-1)  and  (C.3-2)  into  Eq.  (C.4-4),  the  total 
attack  (with  no  trim)  becomes 


a2  » e 


2AjX 

V_ 


+ e 


2A2* 


2A0* 


+ 2KjK2e 


COS 


te--) 


(C.4-5) 


where  t*  has  been  replaced  by  X/V0.  The  trim  angle  of  attack 
is  neglected  so  5*2  may  be  expressed  in  a concise,  workable  form. 
This  is  a reasonable  approximation  since  in  reality  the  trim 
angles  are  small  and  do  not  have  an  appreciable  effect  on  the 
drag. 


Substitution  of  Eq.  (C.4-5)  into  Eq.  (C.4-3),  yields 
the  drag  coefficient 


Kr 


CD  " CD.  + + (CNfl  + «cx) 


2A 


lx 


2A2  x 


Ki2  e + K2e 


2A0x 


+ 2KjK2e  ° cox 


(C.4-6) 


The  downrange  equation  of  motion  is 


« a2x  pw2  y 


* di2 


CD  + ~T  + (CN„  + Kr) 
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2A, 


K2  e V°  + K2  e ° 
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Ji 


(C.4-7 ) 


which  cannot  be  solved  in  closed  form.  As  a result,  perturbations 
are  calculated  for  use  as  corrections.  The  perturbations  are 
based  on  mean  values  of  the  drag  coefficient  with  and  without 
oscillations,  where  the  mean  drag  coef f icier i is  defined  as 


(C.4-8) 


CD 


For  the  oase  without  oscillations  (zero  angle  of  attack) 


The  equation  of  motion  becomes 


(C.4-9) 


« dV  w dV 

g dt  g dx 


or  equivalently 


(C.4-10) 


(C.4-11) 


Substituting  Bq.  (C.4-11)  into  the  Eq.  (C.4-9)  and  integrating 
yields 


Ja  - 0 


1 2 w 

* P A* 


'o 

V 


(C.4-12) 


For  the  case  with  oscillatory  motion,  the  mean  drag 
coefficient  is 


0 


x 


a2  d* 


(C.4-13) 


Upon  substituting  the  expression  for  T2,  Eq.  (C.4-7),  and  inte- 
grating 
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(C.4-14) 
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again  taking  note  of  the  quadrant  for  <t>0  . 

The  mean  drag  coefficients  defined  by  Eqs.  (C.4-12)  and 
(C.4-14)  are  then  used  in  constant  drag  coefficient  equations 
for  the  motion t 


_P  Ag  _ 

" Crs  I 
V - V „e  2 w U 


±L±±(.^  C° * -,) 

pAg  V0  CD  \ ' 


(C.4-16) 


(C.4-17) 


The  two  results  are  subtracted  to  confute  the  perturbation 
quantities  A V and  A t : 


Av  * va  - va-0  - v0 


[ P Ag  ~ pAg  1 

e’  2w  D«X  _ e 2w  Da  - 0 * J 
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For  values  of  range,  X,  prior  to  angle  of  attack  convergence, 
theee  perturbations  are  applied  directly  to  the  velocity  and  flight 
tine  computed  from  the  particle  trajectory  solution  given  in 
Section  C.2,  Eqs.  (C.2-2)  and  (C.2-3).  At  values  of  range  be- 
yond the  point  of  convergence,  the  particle  trajectory  solution 
from  Eqs.  (C.2-2)  and  (C.2-3)  is  restarted  at  the  range  of  con- 
vergence, Xc.  The  initial  conditions  corresponding  to  the  vel- 
ocity and  time  at  tyie  convergence  range,  Xc,  are  the  values  for 
the  particle  trajectory  equations  Without  angle  of  attack  oscil- 
lations) perturbed  by  A V and  A t . 

The  range  at  convergence,  Xc,  is  found  by  determining 
the  point  at  which  the  angle  of  attack  oscillations  decay  to  a 
point  where  the  drag  is  no  longer  appreciably  affected  by  angle 
of  attack.  From  Eq.  (C.4-5),  the  upper  envelope  of  the  oscilla- 
tions is  given  approximately  by 


“UPPER  ' 


(C.4-26) 


The  damping  is  assumed  to  be  small  compared  to  the  natural  fre- 
quency, so  that  many  oscillations  are  required  to  reach  half 
amplitude.  For  the  purposes  of  numerical  calculations,  experi- 
ence has  shown  pitch  oscillations  beyond  the  point  where 
“upper  “ 0.5°  have  little  effect  on  the  trajectory.  Thus,  the 
range  at  convergence  is  defined  as  the  point  at  which  this 
occurs.  (This  is  a computer  code  input  parameter  and  can  be 
easily  changed  to  any  other  value.) 

The  downrange  error,  A X,  at  nominal  time,  tn,  is  com- 
puted by  iterating  the  solution  in  range  until  the  computed 
flight  time  matches  tn.  If  desired,  the  iteration  can  be  avoid- 
ed by  an  additional  approximation.  The  downrange  error  is  close- 
ly approximated  by 


An  ? (tn  - t)  V 


(C.4-21) 


C.4.2  Verification 

Verification  of  the  downrange  perturbation  equations 
was  established  by  comparison  to  four  numerical  solutions  of  the 
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full  aquation*  of  Motion  givan  in  Table  C-3.  The  downrange 
position  agree*  to  within  0.05%.  The  downrange  perturbation, 

A X (the  difference  between  the  downrange  position  with  oscilla 
tory  Motion  and  the  downrange  position  without  oscillatory 
notion  at  the  sssm  flight  tine)  agrees  to  approximately  10%. 


APPENDIX  D 
HITS  PROGRAM  LISTING 


This  appendix  presents  the  HITS  program  listing.  The 
code  was  designed  to  facilitate  conversion  to  real-time 
operation  from  remote  key  board  terminals.  The  program 
is  written  in  FORTRAN  IV  and  requires  less  than  175  K bytes 
of  computer  memory  on  an  IBM  360. 

The  code  generates  seven  (7)  warning  level  (i.e.. 

Level  4)  diagnostics  when  compiled  on  an  IBM- 360/75: 
four  (4)  in  Subroutine  DOABC  and  three  (3)  in  Subroutine 
D0234.  These  are  to  be  expected  and  ignored. 


®oo-«Ni*>*ir>«e  »*•«©>© -tun#  io'Cncoo 
ooooooooo«M«MHMPiai«^NNNNNNMMNNronnnnni<)r)nn't«4  4«<«4««iA 
oeoooooooooooeooooooooooooooooooooaooooooooooooooo 

OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO 

ininininininintnininininintnininwininininininininininininininininininininininininininininintnintntn 

ixxxiixiiiziitiitrxixitiiiitiiiiiizzttiiiittxiiitt 


« 

* 

* 

# 

» 

# 

# 

UJ 

UJ 

* 

» 

« 

UJ 

UJ 

* 

> 

* 

• 

UJ 

UJ 

» 

a 

# 

« 

UJ 

UJ 

UJ 

* 

a 

* 

* 

UJ 

UJ 

UJ 

• 

« 

• 

UJ 

UJ 

UJ 

# 

►»  « 

« 

UJ  UJUJUJUJUJUJ 

OOOQQ 

8 § 

O Q 

8 8 

QOQOQQQ 


OOOOOOQ 
O Q 

o o 

8 8 

o o 


UJ 

z w 
O 3 
M x 
in  u 
►1  u < 
> turn 

S8< 


* 

oooaoao 

* 

m m 

« 

* 

u 

z * 

* 

# 

UJ 

UJ  JZ 

* 

0) 

in  in  in  m 

w u 

# 

-» 

H JO 

* 

in 

in 

in 

u u 

♦ 

< 

inujH 

in 

in 

in 

o u 

* 

a 

>*u> 

* 

in 

in 

in 

o u 

* 

in  oz 

K 

# 

in 

in 

in 

u u 

* 

> 

o 

* 

in 

in 

in 

o u 

* 

h- 

a> 

o z 

u 

« 

m in  m in 

tn 

ooouooo 

* 

X 

U - J 

» 

* 

O 

a 

>Ot-> 

o 

* 

* 

UJ 

<<\J* 

UJ 

# 

¥• 

* 

J 

a 

a 

* 

* 

IL 

< 

< 

♦ 

# 

5 

a 

a 

* 

r 

T 

T 

* 

*-4 

UJ 

UJ 

• 

* 

a 

a 

# 

►* 

* 

> 

a 

a 

* 

H 

# 

h 

« 

# 

M 

* 

# 

u 

# 

* 

o 

* 

M 

»-» 

# 

J 

* 

h"4 

M 

* 

UJ 

» 

**  H* 

W W M 

M 1*4 

* 

> 

# 

* 

a 

* 

M 

* 

UJ 

* 

» 

a 

* 

* 

V 

* 

* 

X 

» 

X XXXXIX 

* 

* 

X 

# 

# 

X 

* 

« 

X 

* 

» 

X 

# 

* 

X 

* 

« 

X 

* 

8 

< - 

J O 

CM 

z - 

< « o 

X N 

O W O' 

UJO  M «* 

i»sjo 

< <z  > 

« • zt-  a 

Onw  j 3 

h>  </>  j a 

o </s  ac  *— • 2 

UJ<«  < 

ax  • *■> 

5S?i  , 

r j j uj 

U_jWin  ►* 


10200 

UJWOO 

k<3kk 


* : s : 

H * * * 

3 # * * 

* #####*##*###*##*####**#* ;<•*#*####***##*#*#■*#*## 

HOOOOOUOOOIJOOOOOOOUOOOOOOOOUOUUUOOUOOOOUUUOUUOVOOOO 


D. 1 MAIN  Program  (continued) 
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D.13  Subroutine  D0234  (continued) 
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D. 28  Subroutine  INCON  (concluded) 
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D.47  Subroutine  SAMPLE  (concluded) 
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D.48  Subroutine  SPRNT  (concluded) 
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